herzo Subroutine

subroutine herzo ( n, x, w )

  !*****************************************************************************80
  !
  !! HERZO computes the zeros the Hermite polynomial Hn(x).
  !
  !  Discussion:
  !
  !    This procedure computes the zeros of Hermite polynomial Ln(x)
  !    in the interval [-1,+1], and the corresponding
  !    weighting coefficients for Gauss-Hermite integration.
  !
  !  Licensing:
  !
  !    This routine is copyrighted by Shanjie Zhang and Jianming Jin.  However, 
  !    they give permission to incorporate this routine into a user program 
  !    provided that the copyright is acknowledged.
  !
  !  Modified:
  !
  !    15 July 2012
  !
  !  Author:
  !
  !    Shanjie Zhang, Jianming Jin
  !
  !  Reference:
  !
  !    Shanjie Zhang, Jianming Jin,
  !    Computation of Special Functions,
  !    Wiley, 1996,
  !    ISBN: 0-471-11963-6,
  !    LC: QA351.C45.
  !
  !  Parameters:
  !
  !    Input, integer ( kind = 4 ) N, the order of the polynomial.
  !
  !    Output, real ( kind = 8 ) X(N), the zeros.
  !
  !    Output, real ( kind = 8 ) W(N), the corresponding weights.
  !
  implicit none

  integer ( kind = 4 ) n

  real ( kind = 8 ) f0
  real ( kind = 8 ) f1
  real ( kind = 8 ) fd
  real ( kind = 8 ) gd
  real ( kind = 8 ) hd
  real ( kind = 8 ) hf
  real ( kind = 8 ) hn
  integer ( kind = 4 ) i
  integer ( kind = 4 ) it
  integer ( kind = 4 ) j
  integer ( kind = 4 ) k
  integer ( kind = 4 ) nr
  real ( kind = 8 ) p
  real ( kind = 8 ) q
  real ( kind = 8 ) r
  real ( kind = 8 ) r1
  real ( kind = 8 ) r2
  real ( kind = 8 ) w(n)
  real ( kind = 8 ) wp
  real ( kind = 8 ) x(n)
  real ( kind = 8 ) x0
  real ( kind = 8 ) z
  real ( kind = 8 ) z0
  real ( kind = 8 ) zl

  hn = 1.0D+00 / n
  zl = -1.1611D+00 + 1.46D+00 * sqrt ( real ( n, kind = 8 ) )

  do nr = 1, n / 2

     if ( nr == 1 ) then
        z = zl
     else
        z = z - hn * ( n / 2 + 1 - nr )
     end if

     it = 0

     do

        it = it + 1
        z0 = z
        f0 = 1.0D+00
        f1 = 2.0D+00 * z
        do k = 2, n
           hf = 2.0D+00 * z * f1 - 2.0D+00 * ( k - 1.0D+00 ) * f0
           hd = 2.0D+00 * k * f1
           f0 = f1
           f1 = hf
        end do

        p = 1.0D+00
        do i = 1, nr - 1
           p = p * ( z - x(i) )
        end do
        fd = hf / p

        q = 0.0D+00
        do i = 1, nr - 1
           wp = 1.0D+00
           do j = 1, nr - 1
              if ( j /= i ) then
                 wp = wp * ( z - x(j) )
              end if
           end do
           q = q + wp
        end do

        gd = ( hd - q * fd ) / p
        z = z - fd / gd

        if ( 40 < it .or. abs ( ( z - z0 ) / z ) <= 1.0D-15 ) then
           exit
        end if

     end do

     x(nr) = z
     x(n+1-nr) = -z
     r = 1.0D+00
     do k = 1, n
        r = 2.0D+00 * r * k
     end do
     w(nr) = 3.544907701811D+00 * r / ( hd * hd )
     w(n+1-nr) = w(nr)

  end do

  if ( n /= 2 * int ( n / 2 ) ) then
     r1 = 1.0D+00
     r2 = 1.0D+00
     do j = 1, n
        r1 = 2.0D+00 * r1 * j
        if ( ( n + 1 ) / 2 <= j ) then
           r2 = r2 * j
        end if
     end do
     w(n/2+1) = 0.88622692545276D+00 * r1 / ( r2 * r2 )
     x(n/2+1) = 0.0D+00
  end if

  return
end subroutine herzo