function r = bisection(fun, xleft, xright, tol, nmax, verbose)

% bisection  Use bisection to find the root of fun(x) = 0
%
% Input:  fun = the function whose zero we are searching
%         xleft, xright = left and right brackets of the root
%         tol = x tolerance
%         nmax = maximum number of iterations
%         verbose = 1 in case all iterations are needed, 0 otherwise
%
% Output: xm = the root of fun(x) = 0 with the desired tolerance

if tol < 3*eps       %  The tolerance is too small
warning(sprintf('the tolerance is too small\n'));
return;
end

if verbose == 1         %  Printing the headings
fprintf(['\niters       left end             midpt' ...
         '             right end        fun(midpt)' ...
	 '  log_2(bracket)\n\n']);
end


%  The main part of the program starts here

a = xleft;            b = xright;
fa = feval(fun,a);    fb = feval(fun,b);

for i = 1:nmax
  deltax = b-a;       xm = a + 0.5*deltax;
  fm = feval(fun,xm);     %  Assigning value fun(xm) to fm

  if verbose == 1         %  Printing all iterations
    fprintf('%3d  %18.14f  %18.14f  %18.14f  %12.3e  %8.3f\n', ...
	     i, a, xm, b, fm, log(a-b)/log(2));
  end

  if deltax < tol
    r = xm;    return;   %  We found the root
  end

  if sign(fm)==sign(fa)
    a = xm;    fa = fm;    %  Root in [xm,b], replace a and fa
  else
    b = xm;    fb = fm;  %  Root in [a,xm], replace b and fb
  end
end

warning(sprintf('root not found within tol after %d iters\n',i));
