ccccccccccccccccccccccccc     Program 3.4     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM S_L_LEGENDRE
C
C Main program for solving the Legendre equation with
C the simplest algorithm for the Sturm-Liouville
C equation and the bisection method for the root search.
C 
      PARAMETER (N=501)
      REAL U(N)
C
C Initialization of the problem
C
      H     =  2.0/(N-1)
      AK    =  0.5
      BK    =  1.5
      DK    =  0.5
      DL    =  1.0E-06
      ISTEP =  0
      U(1)  = -1.0
      U(2)  = -1.0+H
      EK    =  AK
      CALL SMPL (N,H,EK,U)
      F0 = U(N)-1.0
C
C Bisection method for the root
C
      DO      200  WHILE (ABS(DK).GT.DL)
        EK = (AK+BK)/2.0
        CALL SMPL (N,H,EK,U)
        F1 = U(N)-1.0
        IF ((F0*F1).LT.0) THEN
          BK = EK
          DK = BK - AK
        ELSE
          AK = EK
          DK = BK - AK
          F0 = F1
        END IF
        ISTEP = ISTEP + 1
  200 END DO
C
      WRITE (6,999) ISTEP,EK,DK,F1
  999 FORMAT (I4,3F16.8)
      STOP
      END
C
      SUBROUTINE SMPL (N,H,EK,U)
C
C The simplest algorithm for the Sturm-Liouville
C equation.
C
      REAL U(N)
C
      H2 = 2.0*H*H
      Q = EK*(1.0+EK)
      DO      100  I = 2, N-1
        X  =  (I-1)*H-1.0
        P  =  2.0*(1.0-X*X)
        P1 = -2.0*X*H
        U(I+1) = ((2.0*P-H2*Q)*U(I)+(P1-P)*U(I-1))/(P1+P)
  100 CONTINUE
      RETURN
      END