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

%*****************************************************************************80
%
%% regula_falsi2 carries out a modified regula falsi method.
%
%  Discussion:
%
%    Define alpha as the ratio of the new and old interval widths.
%    If the previous step gave alpha > 0.5, take a bisection step.
%    Otherwise, take a regula-falsi (secant) step.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    12 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) < 0, f(xp) > 0.
%
%    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.
%    The test it <= itmax is made.
%
%  Output:
%
%    real xn, xp: real xn, xp: a smaller interval, 
%    for which f(xn) < 0, f(xp) > 0.
%
%    integer it: the number of steps.
%
  it = 0;

  xold = xn;
  x = xp;
  alpha = 0.5;

  while ( true )

    it = it + 1;

    bisect = ( 0.50 < alpha );

    xold = x;
    if ( bisect )
      x = ( xn + xp ) / 2.0;
    else
      x = ( xn * f(xp) - xp * f(xn) ) / ( f(xp) - f(xn) );
    end
%
%  If f(x) < 0.0, update xn.
%
    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
