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

%% bisection2 carries out the bisection method to solve f(x)=0.
%
% 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.
%
  if ( sign ( f(a) ) * sign ( f(b) ) ~= -1 )
    error ( 'F(A) and F(B) are not of opposite signs!' );
  end

  while ( xtol < b - a )
%
%  Get the interval midpoint.
%
    c = ( a + b ) / 2.0;
%
%  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
