% bisection method looks for a zero of function f on the interval
% [a,b].  It stops when the interval length < 2 tol. 
%
% function root = bisect(a,b,tol)

function root = bisect(a,b,tol)

c = (a+b)/2;
i = 0;
while (b-c) > tol
  test = sign(f(b))*sign(f(c));
  if test <= 0 
    a = c;
  else
    b = c;
  end
  c = (a+b)/2;
  i = i+1;
end
i

root = c;