function root = bisection( func ,xl,xu, es, maxit, verbose)

% BISECTION: root = bisection(func,xl,xu,es,maxit):
%   uses bisection method to find the root of a function
% input:
%   func    = name of function 
%   xl, xu  = lower and upper guesses
%   es      = (optional) stopping criterion (%)
%   maxit   = (optional) maximum allowable iterations
%   verbose = (optional) flag for verbose output
% output:
%   root = real root

if func(xl)*func(xu) > 0         % if guesses do not bracket a sign change
   disp('ERROR: No Bracket')     % display an error message
   return                        % and terminate 
end

% if necessary, assign default values

if nargin < 6, verbose=0; end    %if verbose blank set to false...
if nargin < 5, maxit=50; end     %if maxit blank set to 50
if nargin < 4, es=0.001; end     %if es blank set to 0.001

% bisection

iter = 0;
xr   = xl;

% Main loop for bisection algorithm ....

while (1)

  % update solution estimate and iteration no ....

  xrold = xr;
  xr    = (xl + xu)/2;
  iter  = iter + 1;

  if verbose ~= 0,
     fprintf('*** iter = %3d xr = %8.6f\n', iter, xr );
  end;

  % compute relative error

  if xr ~= 0,
     ea = abs((xr - xrold)/xr) * 100;
  end

  % compute update for next point ....

  test = func(xl)*func(xr);
  if test < 0
    xu = xr;
  elseif test > 0
    xl = xr;
  else
    ea = 0;
  end;

  if ea <= es | iter >= maxit,
     break,
  end;
end;

root = xr;