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

%*****************************************************************************80
%
%% regula_falsi3 carries out the modified regula falsi method.
%
%  Discussion:
%
%     The "Illinois algorithm" is used to handle cases where regula falsi
%     seems to be reducing one side of the interval too frequently, at the
%     expense of the other.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    13 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;

  fn = f ( xn );
  fp = f ( xp );
  x = xn;
  fx = f ( xn );

  while ( true )

    it = it + 1;
    xold = x;
    fxold = fx;
    x = ( xn * fp - xp * fn ) / ( fp - fn );
    fx = f ( x );
    old = abs ( xp - xn );
%
%  If f(x) < 0.0, update xn.
%  If we updated xn twice in a row, 
%  shrink fp to encourage updating xp soon.
%
    if ( fx < 0.0 )
      xn = x;
      fn = fx;
      if ( fxold < 0.0 )
        fp = fp / 2.0;
      end
%
%  If f(x) > 0.0, update xp.
%  If we have updated xp twice in a row, 
%  shrink fn to encourage updating xn soon.
%
    else
      xp = x;
      fp = fx;
      if ( fxold > 0.0 )
        fn = fn / 2.0;
      end
    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
