%
% Use bisection to approximate the root of a function.
%
% f(xlo) and f(xhi) must have opposite signs!
%
% The loop terminates when |f(xmid)| < ftol.
% It might also be wise to check if xlo and xhi
% become very close.  That is, stop when |xlo-xhi| < xtol.
%
function x = bisect(f,xlo,xhi,ftol)

    flo  = feval(f,xlo);
    fhi  = feval(f,xhi);
    xmid = (xlo+xhi)/2;
    fmid = feval(f,xmid);
    if (sign(flo*fhi) ~= -1)
        error('bisect: f(xlo) and f(xhi) do not have opposite signs!')
    end
    while (abs(fmid) > ftol)
        if (sign(flo) == sign(fmid))
            xlo = xmid;
            flo = fmid;
        else
            xhi = xmid;
            fhi = fmid;
        end
        xmid = (xlo+xhi)/2;
        fmid = feval(f,xmid);
    end
    x = xmid;