function [ xnew, fxnew ] = zero_muller ( func, fatol, itmax, x1, x2, x3, ...
  xatol, xrtol )

%*****************************************************************************80
%
%% zero_muller() carries out Muller's method, using complex arithmetic.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    26 March 2024
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Gisela Engeln-Muellges, Frank Uhlig,
%    Numerical Algorithms with C,
%    Springer, 1996,
%    ISBN: 3-540-60530-4,
%    LC: QA297.E56213.
%
%  Input:
%
%    fx = FUNC ( x ): evaluates the function.
%
%    real FATOL, the absolute error tolerance for F(X).
%
%    integer ITMAX, the maximum number of steps allowed.
%
%    complex X1, X2, X3, three distinct points to start the iteration.
%
%    real XATOL, XRTOL, absolute and relative
%    error tolerances for the root.
%
%  Output:
%
%    complex XNEW, the estimated root.
%
%    complex FXNEW, the value of the function at XNEW.
%
  xnew = x1;
  xmid = x2;
  xold = x3;

  fxnew = func ( xnew );
  fxmid = func ( xmid );
  fxold = func ( xold );

  fprintf ( 1, '\n' );
  fprintf ( 1, 'zero_muller():\n' );
  fprintf ( 1, '  Muller''s root-finding method (complex root version)\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Iteration     x_real              x_imag             ||fx||           ||disc||\n' );
  fprintf ( 1, '\n' );

  iterate = -2;
  fprintf ( 1, '%6d,%20.10f,%20.10f,%20.10f\n', ...
    iterate, real ( xold ), imag ( xold ), abs ( fxold ) );
  iterate = -1;
  fprintf ( 1, '%6d,%20.10f,%20.10f,%20.10f\n', ...
    iterate, real ( xmid ), imag ( xmid ), abs ( fxmid ) );
  iterate = 0;
  fprintf ( 1, '%6d,%20.10f,%20.10f,%20.10f\n', ...
    iterate, real ( xnew ), imag ( xnew ), abs ( fxnew ) );

  if ( abs ( fxnew ) < fatol )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'zero_muller():\n' );
    fprintf ( 1, '  |F(X)| is below the tolerance.\n' );
    return
  end

  while ( true )
%
%  You may need to swap (XMID,FXMID) and (XNEW,FXNEW).
%
    if ( abs ( fxmid ) <= abs ( fxnew ) )

      c8_temp = xnew;
      xnew = xmid;
      xmid = c8_temp;

      c8_temp = fxnew;
      fxnew = fxmid;
      fxmid = c8_temp;

    end

    xlast = xnew;
    iterate = iterate + 1;

    if ( itmax < iterate )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'zero_muller(): Warning!\n' );
      fprintf ( 1, '  Maximum number of steps taken.\n' );
      break;
    end

    a =  ( ( xmid - xnew ) * ( fxold - fxnew ) ...
         - ( xold - xnew ) * ( fxmid - fxnew ) );

    b = ( ( xold - xnew )^2 * ( fxmid - fxnew ) ...
        - ( xmid - xnew )^2 * ( fxold - fxnew ) );

    c = ( ( xold - xnew ) * ( xmid - xnew ) * ( xold - xmid ) * fxnew );

    xold = xmid;
    xmid = xnew;
%
%  Apply the quadratic formula to get roots XPLUS and XMINUS.
%
    discrm = b * b - 4.0 * a * c;

    if ( a == 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'zero_muller(): Warning!\n' );
      fprintf ( 1, '  The algorithm has broken down.\n' );
      fprintf ( 1, '  The quadratic coefficient A is zero.\n' );
      break;
    end

    xplus = xnew + ( ( - b + sqrt ( discrm ) ) / ( 2.0 * a ) );

    fplus = func ( xplus );

    xminus = xnew + ( ( - b - sqrt ( discrm ) ) / ( 2.0 * a ) );

    fminus = func ( xminus );
%
%  Choose the root with smallest function value.
%
    if ( abs ( fminus ) < abs ( fplus ) )
      xnew = xminus;
    else
      xnew = xplus;
    end

    fxold = fxmid;
    fxmid = fxnew;
    fxnew = func ( xnew );
    fprintf ( 1, '%6d,%20.10f,%20.10f,%20.10f,%20.10f\n', ...
      iterate, real ( xnew ), imag ( xnew ), abs ( fxnew ), abs ( discrm ) );
%
%  Check for convergence.
%
    x_ave = abs ( xnew + xmid + xold ) / 3.0;
    x_inc = xnew - xmid;

    if ( abs ( x_inc ) <= xatol )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'zero_muller():\n' );
      fprintf ( 1, '  Absolute convergence of the X increment.\n' );
      break;
    end

    if ( abs ( x_inc ) <= xrtol * x_ave )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'zero_muller():\n' );
      fprintf ( 1, '  Relative convergence of the X increment.\n' );
      break;
    end

    if ( abs ( fxnew ) <= fatol )
      fprintf ( 1, '\n' );
      fprintf ( 1, 'zero_muller():\n' );
      fprintf ( 1, '  Absolute convergence of |F(X)|.\n' );
      break;
    end

  end

  return
end

