ChebyshevGauss Subroutine

  subroutine ChebyshevGauss(N,nodes,weights)
    implicit none
    integer       :: N
    real(real64)    :: nodes(0:N)
    real(real64)    :: weights(0:N)
    ! Local
    integer    :: j

    do j = 0,N

      weights(j) = pi/(real(N,real64)+1.0_real64)
      nodes(j) = -cos(pi*(2.0_real64*real(j,real64)+1.0_real64)/(2.0_real64*real(N,real64)+2.0_real64))

    enddo

  endsubroutine ChebyshevGauss