!-----  Ex8_9: ルジャンドル多項式   -----
   RECURSIVE FUNCTION P(n,x) RESULT(pnx)
      REAL :: pnx
      INTEGER, INTENT(IN) :: n
      REAL, INTENT(IN) :: x
      IF(n == 0) THEN
         pnx = 1.0
      ELSE IF(n == 1) THEN
         pnx = x
      ELSE
         pnx = (2.0-1.0/n) * x * P(n-1,x) - (1.0-1.0/n) * P(n-2,x)
      END IF
   END FUNCTION
!
   PROGRAM main
      INTEGER :: n, i
      CHARACTER(LEN=1) :: c(-25:25)
      REAL :: x
      PRINT*,'何次の多項式にしますか？'
      READ*,n
      PRINT '( "Legendre Polynomial P(", I2, ", x)" )' , n
      DO i = 0, 20
         x = 0.05*i
         y = p(n, x)
         CALL presentation
         PRINT '(F5.2, 2X, F10.6, 1X, 51A1)', x, y, c
      END DO
   CONTAINS
      SUBROUTINE presentation
         INTEGER :: j
         IF(i /= 0) THEN
             c = ' '
         ELSE
             c = (/ ('-', j = -25, 25) /)
         END IF
         c(-25) = 'I'; c(0) = 'I'; c(25) = 'I'
         c(NINT(25*y)) = '*'
      END SUBROUTINE
   END