function [ x, a, b ] = bisection5 ( xtol, a, b, f )

%% bisection5 carries out the bisection method to solve f(x)=0.
%
% This code will stop if too many bisections are tried.
%
% xtol, tolerance for interval size.
% a, b, the endpoints of a change of sign interval.
%
% x, output, the estimated solution.
% a, b, output, the updated endpoints of the change of sign interval.
%
  global bisect_num;

  fa = sign ( f(a) );

  bisect_max = 50;
  bisect_num = 0;

  while ( xtol < b - a )
%
%  Get the interval midpoint.
%
    c = ( a + b ) / 2.0;
    bisect_num = bisect_num + 1;

    if ( bisect_max < bisect_num )
      fprintf ( 'Maximum number of steps exceeded!\n' );
      fprintf ( 'Bisection terminated without satisfying tolerance.\n' );
      x = ( a + b ) / 2.0;
      return
    end
%
%  Choose new interval [C,B] or [A,C].
%
    if ( sign ( f(c) ) == sign ( fa ) )
      a = c;
    else
      b = c;
    end

  end
%
%  Choose midpoint as estimated zero.
%
  x = ( a + b ) / 2.0;

  return
end
