function root = bisect(fname, interval, tol);

% BISECT  Bisection method
%
% Input       fname         : Name of function.
%             interval      : Interval containing the root. 
%             tol           : Absolute error of the root.
% Output      root          : Position of the root.
%
% Remarks  1) The root must be in the interval (lower, upper).
%          2) If the function has more roots in this interval,
%             it will not be clear which root will be found
%          3) As the initial interval is halved at each iteration, it is easy 
%             to calculate the necessary number of iterations to match the
%             required tolerance.

if nargin < 3, tol = 1e-7; end

lower = interval(1);
upper = interval(2);
niter = ceil(log((upper-lower)/tol)/log(2));

for n=1:niter,
  root = (lower+upper)/2;

% Determine which half of the interval contains the root and half the 
% interval

  flower = feval(fname,lower);
  froot = feval(fname,root);
  if flower*froot > 0,
    lower = root;
  else
    upper = root;
  end
end