subroutine zero_muller ( func, fatol, itmax, x1, x2, x3, xatol, 
     &  xrtol, xnew, fxnew )

c*********************************************************************72
c
cc zero_muller() carries out Muller's method, using complex arithmetic.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    27 March 2024
c
c  Author:
c
c    John Burkardt
c
c  Reference:
c
c    Gisela Engeln-Muellges, Frank Uhlig,
c    Numerical Algorithms with C,
c    Springer, 1996,
c    ISBN: 3-540-60530-4,
c    LC: QA297.E56213.
c
c  Input:
c
c    external FUNC, the name of the routine that evaluates the function.
c    FUNC should have the form:
c      subroutine func ( x, fx )
c      double complex fx
c      double complex x
c
c    double precision FATOL, the absolute error tolerance for F(X).
c
c    integer ITMAX, the maximum number of steps allowed.
c
c    double complex X1, X2, X3, three distinct points to start the
c    iteration.
c
c    double precision XATOL, XRTOL, absolute and relative
c    error tolerances for the root.
c
c  Output:
c
c    double complex XNEW, the estimated root.
c
c    double complex FXNEW, the value of the function at XNEW.
c
      implicit none

      double complex a
      double complex b
      double complex c
      double complex c8_temp
      double complex discrm
      double precision fatol
      double complex fminus
      double complex fplus
      external func
      double complex fxmid
      double complex fxnew
      double complex fxold
      integer iterate
      integer itmax
      double precision x_ave
      double complex x_inc
      double complex x1
      double complex x2
      double complex x3
      double precision xatol
      double complex xlast
      double complex xmid
      double complex xminus
      double complex xnew
      double complex xold
      double complex xplus
      double precision xrtol

      xnew = x1
      xmid = x2
      xold = x3

      call func ( xnew, fxnew )
      call func ( xmid, fxmid )
      call func ( xold, fxold )

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'zero_muller():'
      write ( *, '(a)' ) '  Muller''s root-finding method.'
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  Iteration     x_real              ' // 
     &  'x_imag             ||fx||           ||disc||'
      write ( *, '(a)' ) ' '

      iterate = -2
      write ( *, '(i6,f20.10,f20.10,f20.10)' ) 
     &  iterate, xold, abs ( fxold )
      iterate = -1
      write ( *, '(i6,f20.10,f20.10,f20.10)' ) 
     &  iterate, xmid, abs ( fxmid )
      iterate = 0
      write ( *, '(i6,f20.10,f20.10,f20.10)' ) 
     &  iterate, xnew, abs ( fxnew )

      if ( abs ( fxnew ) < fatol ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'zero_muller():'
        write ( *, '(a)' ) '  |F(X)| is below the tolerance.'
        return
      end if

      do
c
c  You may need to swap (XMID,FXMID) and (XNEW,FXNEW).
c
        if ( abs ( fxmid ) <= abs ( fxnew ) ) then

          c8_temp = xnew
          xnew = xmid
          xmid = c8_temp

          c8_temp = fxnew
          fxnew = fxmid
          fxmid = c8_temp

        end if

        xlast = xnew
        iterate = iterate + 1

        if ( itmax < iterate ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'zero_muller(): Warningc'
          write ( *, '(a)' ) '  Maximum number of steps taken.'
          exit
        end if

        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
c
c  Apply the quadratic formula to get roots XPLUS and XMINUS.
c
        discrm = b**2 - 4.0D+00 * a * c

        if ( a == 0.0D+00 ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'zero_muller(): Warningc'
          write ( *, '(a)' ) '  The algorithm has broken down.'
          write ( *, '(a)' ) '  The quadratic coefficient A is zero.'
          exit
        end if

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

        call func ( xplus, fplus )

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

        call func ( xminus, fminus )
c
c  Choose the root with smallest function value.
c
        if ( abs ( fminus ) < abs ( fplus ) ) then
          xnew = xminus
        else
          xnew = xplus
        end if

        fxold = fxmid
        fxmid = fxnew
        call func ( xnew, fxnew )
        write ( *, '(i6,f20.10,f20.10,f20.10,f20.10)' ) 
     &    iterate, xnew, abs ( fxnew ), abs ( discrm )
c
c  Check for convergence.
c
        x_ave = abs ( xnew + xmid + xold ) / 3.0D+00
        x_inc = xnew - xmid

        if ( abs ( x_inc ) <= xatol ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'zero_muller():'
          write ( *, '(a)' ) '  Absolute convergence of X increment.'
          exit
        end if

        if ( abs ( x_inc ) <= xrtol * x_ave ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'zero_muller():'
          write ( *, '(a)' ) '  Relative convergence of X increment.'
          exit
        end if

        if ( abs ( fxnew ) <= fatol ) then
          write ( *, '(a)' ) ' '
          write ( *, '(a)' ) 'zero_muller():'
          write ( *, '(a)' ) '  Absolute convergence of |F(X)|.'
          exit
        end if

      end do

      return
      end