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

%*****************************************************************************80
%
%% whittaker2() implements the double convex acceleeration of Whittaker's method.
%
%  Discussion:
%
%    x(n+1) = x(n) - (1/4) * ( f(x(n)) / f'(x(n)) ) *
%      ( 2 - L(x(n)) + ( 2 + L(x(n)) ) / ( 1 - L(x(n)) + 0.5 * L(x(n))**2 ) )
%
%    where
%
%      L(x) = f(x) * f''(x) / f'(x)**2
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 May 2024
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X: an estimate for the root of the equation.
%
%    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: an approximate root for which abs ( F(X) ) <= ABSERR.
%
%    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.
%
  ierror = 0;
  k = 0;
  fx = f ( x, 0 );
  dfx = f ( x, 1 );
  d2fx = f ( x, 2 );
%
%  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

    if ( dfx == 0.0 )
      ierror = 3;
      return
    end

    u = d2fx * fx / dfx / dfx;

    if ( 1.0 - u + 0.5 * u * u  == 0.0 )
      ierror = 4;
      return
    end
%
%  Set the increment.
%
    dx = - 0.25 * ( fx / dfx ) ...
      * ( 2.0 - u + ( 2.0 + u ) / ( 1.0 - u + 0.5 * u * u ) );
%
%  Update the iterate and function values.
%
    x = x + dx;
    fx = f ( x, 0 );
    dfx = f ( x, 1 );
    d2fx = f ( x, 2 );

  end

  return
end
