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

%*****************************************************************************80
%
%% r8_muller() implements Muller's method for real arithmetic.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    19 July 2007
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, X1, X2: three 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, X2: X is an approximate root satisfying
%    abs ( F(X) ) < ABSERR, and X1 and X2 are the previous estimates.
%
%    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;
  fx2 = f ( x2, 0 );
  fx1 = f ( x1, 0 );
  fx = f ( x, 0 );
%
%  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

    q = ( x - x1 ) / ( x1 - x2 );

    a = q * fx - q * ( 1.0 + q ) * fx1 + q * q * fx2;
    b = ( 2.0 * q + 1.0 ) * fx - ( 1.0 + q ) * ( 1.0 + q ) * fx1 + q * q * fx2;
    c = ( 1.0 + q ) * fx;

    term = b * b - 4.0 * a * c;
    term = max ( term, 0.0 );
    term = sqrt ( term );
    if ( b < 0.0 )
      term = - term;
    end
%
%  Set the increment.
%
    dx = - ( x - x1 ) * 2.0 * c / ( b + term );
%
%  Remember current data for next step.
%
    x2 = x1;
    fx2 = fx1;
    x1 = x;
    fx1 = fx;
%
%  Update the iterate and function values.
%
    x = x + dx;
    fx = f ( x, 0 );

  end

  return
end
