function [ xn, xp, it ] = bisection3 ( f, xn, xp, xtol, ftol, itmax )

%*****************************************************************************80
%
%% bisection3 carries out the bisection method, and computes alpha.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    09 September 2019
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    function handle F: a pointer to a function.  This can be
%    an anonymous function or the name of a MATLAB M-file function.
%
%    real xn, xp: f(xn) is negative, f(xp) is positive.
%
%    real xtol: the interval size tolerance.  
%    The test |xn-xp| <= xtol is made.
%
%    real ftol: the residual tolerance.  
%    The test |f(x)| <= ftol is made.
%
%    integer itmax: the maximum number of steps allowed.
%
%  Output:
%
%    real xn, xp: f(xn) is negative, f(xp) is positive, and
%    the estimated root is the midpoint of this interval.
%
%    integer it: the number of bisection steps.
%
  it = 0;

  while ( true )

    it = it + 1;

    x = ( xn + xp ) / 2.0;

    old = abs ( xp - xn );

    if ( f ( x ) < 0.0 )
      xn = x;
    else
      xp = x;
    end

    new = abs ( xp - xn );
    alpha = new / old;
    fprintf ( 1, 'Step %d: alpha = %g\n', it, alpha );
%
%  Satisfied?
%
    if ( ...
      ( new <= xtol ) && ...
      ( abs ( f ( x ) ) <= ftol ) ...
    )
      return
    end
%
%  Dissatisfied?
%
    if ( itmax <= it )
      return
    end
%
%  Repeat.
%
  end

  return
end
