! ルジャンドル多項式
RECURSIVE FUNCTION p( n, x ) RESULT( pnx )     ! 必ず RESULT を使うこと

  IMPLICIT NONE
  INTEGER, INTENT(IN) :: n
  REAL, INTENT(IN) :: x
  REAL :: pnx

  IF ( n == 0 ) THEN
    pnx = 1.0E0
  ELSE IF ( n == 1 ) THEN
    pnx = x
  ELSE
    pnx = ( 2.0E0 - 1.0E0 / n ) * p(n-1,x) * x &
        - ( 1.0E0 - 1.0E0 / n ) * p(n-2,x)
  END IF

END FUNCTION

! main
PROGRAM main

  IMPLICIT NONE
  INTEGER :: n, i
  REAL :: x, y, p
  CHARACTER(LEN=1) :: c(-25:25)

  PRINT '(1X, A)', '何次の多項式にしますか？'
  READ *, n

  PRINT '("Legendre polynomial p(",I2, ", x)")', n
  DO i = 0, 20
    x = 0.05E0 * i
    y = p(n,x)
    CALL presentation()
    PRINT '(F5.2, 2X, F10.6, 1X, 51A1)', x, y, c
  END DO

CONTAINS

  SUBROUTINE presentation()    !  内部手続き故，変数を運ぶ必要なし

    IMPLICIT NONE
    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 presentation

END PROGRAM main