function r = bisect(f,int,tol)
%
% BISECT(f,a,b,tol)
% bisection method to find a root of f
%
% input:
%   f     string giving function (e.g., 'sin(x^3)+3*x-5')
%   int   vector of length two specifying a bracketing interval
%   tol   desired size of final bracketing interval
%         If tol is too small (e.g., tol=0), the smallest possible
%         interval (based on machine precision) will be found.
%         If tol is omitted, tol=0 is used.
%
%  output:
%    r   approximate root
%           r = (a'+b')/2 where f(a')*f(b') < 0 and abs(b' - a') < tol
%
% This is an implementation of the bisection root-finding algorithm.
% Note: there are much better root-finding algorithms available.
%
% Douglas N. Arnold, September 1997

% check for 2 or 3 arguments, set tol=0 if not given
if (nargin == 2)
  tol = 0;
elseif (nargin ~= 3)
  fprintf('\nUsage: r = root(f,int,tol)\n');
  return
end

% interpret input function so it can be evaluated with feval
f=fcnchk(f);
% initialize a and b and evaluate function there, start with c=a
a=int(1); b=int(2); c=a;
fa = feval(f,a); fb = feval(f,b); fc = fa;

% check for exact root in input data
if (fa == 0)
  r = a; fprintf('\nexact root found:  %18.16g\n',r); return
end
if (fb == 0)
  r = b; fprintf('\nexact root found:  %18.16g\n',r); return
end

% make sure input interval brackets
if (sign(fa) == sign(fb))
  fprintf('\ninput interval does not bracket\n'); return
end

while 1 %%%%% START OF MAIN LOOP
  % check for convergence
    if (abs(a-b) <= tol)
      r = b;
      fprintf('\ntolerance achieved\n');
      return
    end

  % compute midpoint and evaluate function there
    c = (a +  b)/2; fc = feval(f,c);

  % check if midpoint equals endpoint to machine precision;
  % if so stop, because we can't get a smaller interval
    if (c == a | c == b)
      r = c;
      if (tol > 0); fprintf('\ntolerance not achievable\n'); end
      return
    end
    
  % replace bracketing interval with bracketing half interval
    if ( sign(fa) == sign(fc) )
      a = c;
      fa = fc;
    else
      b = c;
      fb = fc;
    end
end %%%%% END OF MAIN LOOP
