!!!!!!!!!!!!!!!!!!!!!!!!!!!   Program 3.4   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                       !
! Please Note:                                                          !
!                                                                       !
! (1) This computer program is written by Tao Pang in conjunction with  !
!     his book, "An Introduction to Computational Physics," published   !
!     by Cambridge University Press in 1997.                            !
!                                                                       !
! (2) No warranties, express or implied, are made for this program.     !
!                                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
PROGRAM S_L_LEGENDRE
!
! Main program for solving the Legendre equation with the simplest 
! algorithm for the Sturm-Liouville equation and the bisection method
! for the root search.  Copyright (c) Tao Pang 1997.
! 
  IMPLICIT NONE
  INTEGER, PARAMETER :: N=501
  INTEGER :: ISTEP
  REAL :: DL,H,AK,BK,DK,EK,F0,F1
  REAL, DIMENSION (N) :: U
!
! Initialization of the problem
!
  DL = 1.0E-06
  H  = 2.0/(N-1)
  AK = 0.5
  BK = 1.5
  DK = 0.5
  EK = AK
  U(1)  = -1.0
  U(2)  = -1.0+H
  ISTEP =  0
  CALL SMPL (N,H,EK,U)
  F0 = U(N)-1.0
!
! Bisection method for the root
!
  DO 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
  END DO
!
  WRITE (6,"(I4,3F16.8)") ISTEP,EK,DK,F1
END PROGRAM S_L_LEGENDRE
!
SUBROUTINE SMPL (N,H,EK,U)
!
! The simplest algorithm for the Sturm-Liouville equation.
! Copyright (c) Tao Pang 1997.
!
  IMPLICIT NONE
  INTEGER, INTENT (IN) :: N
  INTEGER :: I
  REAL, INTENT (IN) :: H,EK
  REAL :: H2,Q,X,P,P1
  REAL, INTENT (OUT), DIMENSION (N) :: U
!
  H2 = 2.0*H*H
  Q = EK*(1.0+EK)
  DO 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)
  END DO
END SUBROUTINE SMPL