function [b, steps] = bisect (fun,x,tol)
% BISECT:  zero of a function of one
% variable via the bisection method.
% bisect (fun,x) returns a zero of the
% function fun.  fun is a function
% handle or a string with the name of a
% function.  x is a starting guess. The
% value of b returned is near a point
% where fun changes sign.  For example,
% bisect (@sin,3) is pi.  Note the use
% of the function handle, @sin.
%
% An optional third input argument sets
% a tolerance for the relative accuracy
% of the result.  The default is eps.
% An optional second output argument
% gives a matrix containing a trace of
% the steps; the rows are of the form
% [c (f(c))].
%
% MATLAB Primer, 6th Edition
% Kermit Sigmon and Timothy A. Davis
% Section 8.1, page 39.

if (nargin < 3)
    % default tolerance
    tol = eps ;
end
trace = (nargout == 2) ;
if (x ~= 0)
    dx = x/20 ;
else
    dx = 1/20 ;
end
a = x - dx ;
fa = feval (fun, a) ;
b = x + dx ;
fb = feval (fun, b) ;
if (trace)
    steps = [a fa ; b fb] ;
end

% find a change of sign
while (fa > 0) == (fb > 0)
    dx = 2*dx ;
    a = x - dx ;
    fa = feval (fun, a) ;
    if (trace)
        steps = [steps ; [a fa]] ;
    end
    if (fa > 0) ~= (fb > 0)
        break
    end
    b = x + dx ;
    fb = feval (fun, b) ;
    if (trace)
        steps = [steps ; [b fb]] ;
    end
end

% main loop
while (abs (b-a) > 2*tol*max(abs(b),1))
    c = a + (b-a)/2 ;
    fc = feval (fun, c) ;
    if (trace)
        steps = [steps ; [c fc]] ;
    end
    if (fb > 0) == (fc > 0)
        b = c ;
        fb = fc ;
    else
        a = c ;
        fa = fc ;
    end
end