function root = bisect(f,a,b,tol)
%
% BISECT(f,a,b,tol)
% bisection method to find a root of f
%
% input:
%   f       string giving function name
%   a,b     such that f(a)*f(b) < 0
%   tol     positive tolerance
%
%  output:
%    root   approximate root
%           root = (a'+b')/2 where f(a')*f(b') < 0 and abs(b' - a') < tol
%
% The function f should be specified in a file f.m of the
% form "function y=f(x)".

% check input data
fa = feval(f,a);
fb = feval(f,b);
if ( fa*fb > 0 )
   fprintf('input interval does not bracket\n');
   return
end
if ( tol <= 0 )
   fprintf('tolerance must be positive\n');
   return
end
% main loop
while 1
% uncomment next line to get printed output each iteration
  fprintf('%18.14f %18.14f %18.14f %18.14f\n',a,b,fa,fb);
% check if tolerance is met	
  if ( abs(b-a) <= tol )
    root = (a+b)/2;
    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)
    fprintf('maximum possible precision achieved\n');
    root = c;
    return
  end
% replace bracketing interval with bracketing half interval
  if ( fa*fc > 0)
    a = c;
    fa = fc;
  else
    b = c;
    fb = fc;
  end
end
