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

% BISECT_DEMO  Bisection method (demo)
%
% 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-3; end
if nargin < 2, interval = [0,1]; end
if nargin < 1, fname = 'fff'; end

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

x = lower:(upper-lower)/99:upper;
y = feval(fname,x);
clf
plot(x,y,x,zeros(size(x)),'--');
title('Bisection Demo (<enter> in command window gives next point)')
h = [];
hold on

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 length(h),
    set(h,'Marker','o','MarkerSize',4)
  end
  h = line(root,froot,'Marker','*','MarkerSize',10,'Color','r');
  pause
  if flower*froot > 0,
    lower = root;
  else
    upper = root;
  end
end
title('Bisection Demo (end)')