function approx_root = bisect ( f, xa, xb )
% bisect(f,xa,xb) computes the root of a function f using bisection
% on the interval [xa,xb].  It stops when the width of the interval
% is reduced below 0.000001

% original written by John Burkardt
% modified M. M. Sussman, 9/16/00
% see Atkinson, Section 2.1

fa = feval ( f, xa ); 
fb = feval ( f, xb );

% A little error checking never hurts

if ( sign ( fa ) == sign ( fb ) )
  'BISECT - Fatal error!'
  '  [A,B] is not a change-of-sign interval.'
  root = 0;
  return
end

if (xb < xa)
  'BISECT - Fatal error!'
  '  left and right endpoints interchanged.'
  root=0;
  return
end

% Now start the real work

while ( abs ( xb - xa ) > 0.000001 )

  xc = ( xa + xb ) / 2;
  approx_root = xc;
  fc = feval ( f, xc );

  % print out iteration values to see how it is going
  fprintf('\nxa=%f,\txc=%f,\txb=%f\n',xa,xc,xb);
  fprintf('fa=%f,\tfc=%f,\tfb=%f\n',fa,fc,fb);

  if ( sign(fb) * sign(fc) <= 0 )
    xa = xc;
    fa = fc;
  else
    xb = xc;
    fb = fc;
  end

end