function [p,iter]=bisection(f,a,b,tol,Nmax)
% Technically this will always converge so we shouldn't need an Nmax
% Ex:
% f = @(x) x.^2 - 5;
% [p,iter] = bisection(f,0,5,1e-6,100)


if nargin == 3
   tol = 1e-6; Nmax = 10^6; 
end
if nargin ==4
    Nmax = 10^6;
end

iter = 1;
while iter < Nmax
  p=(a+b)/2;
  if p-a<tol || abs(f(p)) < 100*eps
      break; 
  end
  if f(a)*f(p)>0
    a=p;
  else
    b=p;
  end
  iter = iter+1;
end

if iter > Nmax
   disp('Maximum Iterations reached...') 
end

end
