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

%% bisection3 carries out the bisection method to solve f(x)=0.
%
% This code will stop if more than 52 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.
% bisect_num, the number of bisection steps.
%
  if ( sign ( f(a) ) * sign ( f(b) ) ~= -1 )
    error ( 'F(A) and F(B) are not of opposite signs!' );
  end

  bisect_max = 52;
  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 ( f(a) ) )
      a = c;
    else
      b = c;
    end

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

  return
end
