function root = bisect(func,x1,x2,filter,tol)
% Finds a bracketed zero of f(x) by bisection.
% USAGE: root = bisect(func,x1,x2,filter,tol)
% INPUT:
% func   = handle of function that returns f(x).
% x1,x2  = limits on interval containing the root.
% filter = singularity filter: 0 = off (default), 1 = on.
% tol    = error tolerance (default is 1.0e4*eps).
% OUTPUT:
% root   = zero of f(x), or NaN if singularity suspected.

if nargin < 5; tol = 1.0e4*eps; end
if nargin < 4; filter = 0; end
f1 = feval(func,x1);
if f1 == 0.0; root = x1; return; end
f2 = feval(func,x2);
if f2 == 0.0; root = x2; return; end
if sign(f1)== sign(f2)
    error('Root is not bracketed in (x1,x2)')
end
n = ceil(log(abs(x2 - x1)/tol)/log(2.0));
for i = 1:n
    x3 = 0.5*(x1 + x2);
    f3 = feval(func,x3);
    if(filter == 1) & (abs(f3) > abs(f1))...
                    & (abs(f3) > abs(f2))
        root = NaN; return
    end
    if f3 == 0.0
        root = x3; return
    end
    if sign(f2) ~= sign(f3)
        x1 = x3; f1 = f3;
    else
        x2 = x3; f2 = f3;
    end
end
root = (x1 + x2)/2;
