!-----  Ex8_9: ルジャンドル多項式   -----
   REAL FUNCTION P(n,x)
      REAL :: pnx(0:n) 
      INTEGER, INTENT(IN) :: n
      REAL, INTENT(IN) :: x
      INTEGER :: i
      pnx(0) = 1.0 ;  pnx(1) = x
      DO i = 2, n
         pnx(i) = (2.0-1.0/i)*x*pnx(i-1) - (1.0-1.0/i)*pnx(i-2)
      END DO
      P = pnx(n)     ! 左辺には，引数の付かない関数名だけを書く
   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