function root = bisect0(f,a,b,tol)
%
% bisection method to find a root of f
% see bisect.m for more complete implementation

fa = feval(f,a); fb = feval(f,b);
while 1
% 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);
% replace bracketing interval with bracketing half interval
  if ( fa*fc > 0)
    a = c;
    fa = fc;
  else
    b = c;
    fb = fc;
  end
end
