function [zero, res, niter] = bisect(fun, a, b, tol, nmax, varargin)
x = [a,(a+b)/2,b];
fx = feval(fun,x,varargin{:});
if fx(1)*fx(3)>0
error('need opposite signs at extremal points')
elseif fx(1) == 0
zero = a;
niter = 0; res = 0;
return
elseif fx(3) == 0
zero = b;
niter = 0; res = 0;
return
end
I = (b-a)/2;
niter = 0;
while I>= tol & niter <= nmax
niter = niter + 1;
if sign(fx(1))*sign(fx(2)) < 0
x(3) = x(2); x(2) = x(1) + (x(3) - x(1))/2;
fx = feval(fun,x,varargin{:}); I = (x(3) - x(1))/2;
elseif sign(fx(2))*sign(fx(3)) < 0
x(1) = x(2); x(2) = x(1) + (x(3) - x(1))/2;
fx = feval(fun,x,varargin{:}); I = (x(3) - x(1))/2;
else
x(2) = x(find(fx==0)); I = 0;
end
end
if niter > nmax
fprintf(['maximum iterations reached without achieving desired tolerance'])
end
zero = x(2); x = x(2); res = feval(fun,x,varargin{:});