function [ x, x1, ierror, k ] = rhein2 ( x, x1, abserr, kmax, f )

%*****************************************************************************80
%
%% rhein2() implements the Rheinboldt bisection/secant/inverse quadratic method.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Werner Rheinboldt,
%    Algorithms for Finding Zeros of a Function,
%    UMAP Journal,
%    Volume 2, Number 1, 1981, pages 43 - 72.
%
%  Input:
%
%    real X, X1: two distinct points that start the method.
%
%    real ABSERR, an error tolerance.
%
%    integer KMAX, the maximum number of iterations allowed.
%
%    real external F, the name of the routine that
%    evaluates the function or its derivatives, of the form
%      function f ( x, ider )
%
%  Output:
%
%    real X, X1: X is an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 is the previous estimate.
%
%    integer IERROR, error indicator.
%    0, no error occurred.
%    nonzero, an error occurred, and the iteration was halted.
%
%    integer K, the number of steps taken.
%

%
%  Initialization.
%
  i = 0;
  ierror = 0;
  k = 0;
  fx = f ( x, 0 );
  fx1 = f ( x1, 0 );

  if ( abs ( fx1 ) < abs ( fx ) )
    s = x;
    x = x1;
    x1 = s;
    s = fx;
    fx = fx1;
    fx1 = s;
  end

  x2 = x1;
  fx2 = fx1;
  k = 0;
  t = 0.5 * abs ( x - x1 );
%
%  Iteration loop:
%
  while ( true )
%
%  If the error tolerance is satisfied, then exit.
%
    if ( abs ( fx ) <= abserr )
      break
    end

    k = k + 1;

    if ( kmax < k )
      ierror = 2;
      return
    end

    i = i + 1;

    em = 0.5 * ( x1 - x );
%
%  Compute the numerator and denominator for secant step.
%
    if ( 2.0 * abs ( x2 - x ) < abs ( x1 - x ) )
      ps = ( x - x2 ) * fx;
      qs = fx2 - fx;
    else
      ps = ( x - x1 ) * fx;
      qs = fx1 - fx;
    end

    if ( ps < 0.0 )
      ps = - ps;
      qs = - qs;
    end
%
%  Compute the numerator and denominator for inverse quadratic.
%
    piq = 0.0;
    qiq = 0.0;

    if ( x1 ~= x2 )
      u = fx / fx2;
      v = fx2 / fx1;
      w = fx / fx1;
      piq = u * ( 2.0 * em * v * ( v - w ) - ( x - x2 ) * ( w - 1.0 ) );
      qiq = ( u - 1.0 ) * ( v - 1.0 ) * ( w - 1.0 );

      if ( 0.0 < piq )
        qiq =  - qiq;
      end

      piq = abs ( piq );
    end
%
%  Save the old minimum residual point.
%
    x2 = x;
    fx2 = fx;

    stpmin = ( abs ( x ) + abs ( em ) + 1.0 ) * abserr;
%
%  Choose bisection, secant or inverse quadratic step.
%
    forced = false;

    if ( 3 < i )
      if ( t < 8.0 * abs ( em ) )
        forced = true;
      else
        i = 0;
        t = em;
      end
    end
%
%  Set the increment.
%
    if ( forced )
      dx = em;
    elseif ( piq < 1.5 * em * qiq && abs ( qiq ) * stpmin < abs ( piq ) )
      dx = piq / qiq;
    elseif ( ps < qs * em && abs ( qs ) * stpmin < abs ( ps ) )
      dx = ps / qs;
    else
      dx = em;
    end
%
%  Update the iterate and function values.
%
    x = x + dx;
    fx = f ( x, 0 );
%
%  Set the new X1 as either X1 or X2, depending on whether
%  F(X1) or F(X2) has the opposite sign from F(X).
%
    if ( sign ( fx ) == sign ( fx1 ) )
      x1 = x2;
      fx1 = fx2;
    end
%
%  Force ABS ( FX ) <= ABS ( FX1 ).
%
    if ( abs ( fx1 ) < abs ( fx ) )
      s = x;
      x = x1;
      x1 = s;
      s = fx;
      fx = fx1;
      fx1 = s;
    end

  end

  return
end
