ccccccccccccccccccccccccc     Program 3.D     ccccccccccccccccccccccccccc
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                       c
c Please Note:                                                          c
c                                                                       c
c (1) This computer program is written by Tao Pang in conjunction with  c
c     his book, "An Introduction to Computational Physics," published   c
c     by Cambridge University Press in 1997.                            c
c                                                                       c
c (2) No warranties, express or implied, are made for this program.     c
c                                                                       c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
                    PROGRAM SCHR
C
C Main program for solving the eigenvalue problem of the
C one-dimensional Schroedinger equation. 
C Copyright (c) Tao Pang 1997.
C
      PARAMETER (N=501,M=5,IMAX=100)
      REAL UL(N),UR(N),QL(N),QR(N),S(N)
C
      H2M   =  0.5 
      EA    =  2.4 
      EB    =  2.7
      E     =  EA
      DE    =  0.1
      XL0   = -10.0
      XR0   =  10.0
      H     =  (XR0-XL0)/(N-1)
      C     =  1.0/H2M
      DL    =  1.0E-06
      UL(1) =  0.0
      UL(2) =  0.01
      UR(1) =  0.0
      UR(2) =  0.01
C
C Set up the potential q(x) and source s(x)
C
      DO    100   I = 1, N
        XL    = XL0+(I-1)*H
        XR    = XR0-(I-1)*H
        QL(I) = C*(E-V(XL))
        QR(I) = C*(E-V(XR))
        S(I)  = 0.0
  100 CONTINUE 
C
C Find the matching point at the right turning point
C
      DO     200  I = 1, N-1
        IF(((QL(I)*QL(I+1)).LE.0).AND.(QL(I).GT.0)) IM = I
  200 CONTINUE 
C    
C Numerov algorithm from left to right and vice versa
C
      NL = IM+1
      NR = N-IM+2
      CALL NMRV2 (NL,H,QL,S,UL)
      CALL NMRV2 (NR,H,QR,S,UR)
C
C Rescale the left solution
C
      FACT = UR(NR-1)/UL(IM)
      DO      300  I = 1, NL
        UL(I) = FACT*UL(I)
  300 CONTINUE 
C
      F0 = UR(NR)+UL(NL)-UR(NR-2)-UL(NL-2)
      F0 = F0/(2.0*H*UR(NR-1))
C
C Bisection method for the root
C
      ISTEP = 0
      DO      400  WHILE ((ABS(DE).GT.DL).AND.(ISTEP.LT.IMAX))
        E1 = E
        E  = (EA+EB)/2.0
        DO    310   I = 1, N
          QL(I) = QL(I)+C*(E-E1)
          QR(I) = QR(I)+C*(E-E1)
  310   CONTINUE 
C
C Find the matching point at the right turning point
C
        DO   320  I = 1, N-1
          IF(((QL(I)*QL(I+1)).LE.0).AND.(QL(I).GT.0)) IM = I
  320   CONTINUE 
C    
C Numerov algorithm from left to right and vice versa
C
        NL = IM+1
        NR = N-IM+2
        CALL NMRV2 (NL,H,QL,S,UL)
        CALL NMRV2 (NR,H,QR,S,UR)
C
C Rescale the left solution
C
        FACT = UR(NR-1)/UL(IM)
        DO    330  I = 1, NL
          UL(I) = FACT*UL(I)
  330   CONTINUE 
C
        F1 = UR(NR)+UL(NL)-UR(NR-2)-UL(NL-2)
        F1 = F1/(2.0*H*UR(NR-1))
C
        IF ((F0*F1).LT.0) THEN
          EB = E 
          DE = EB-EA
        ELSE
          EA = E
          DE = EB-EA
          F0 = F1
        END IF
        ISTEP = ISTEP+1
  400 END DO
C
      SUM = 0.0
      DO    450  I = 1, N
        IF(I.GT.IM) UL(I) = UR(N-I)
        SUM = SUM+UL(I)*UL(I)
  450 CONTINUE
      SUM=SQRT(H*SUM)

      WRITE(6,9998) ISTEP,IMAX
      WRITE(6,9999) E,DE,F0,F1
      DO    500 I = 1, N, M
        XL = XL0+(I-1)*H
        UL(I) = UL(I)/SUM
        WRITE(15,9999) XL,UL(I)
        WRITE(16,9999) XL,V(XL)
  500 CONTINUE
C
      STOP
 9998 FORMAT(2I4)
 9999 FORMAT(4F20.8)
      END
C
      FUNCTION V(X)
        A = 1.0
        B = 4.0
        V = 3.0-A*A*B*(B-1.0)/(COSH(A*X)**2)/2.0
      RETURN
      END
C
      SUBROUTINE NMRV2 (N,H,Q,S,U)
C
C The Numerov algorithm for the equation u"(x)+q(x)u(x)=s(x)
C as given in Eqs. (3.82)-(3.85) in the book.
C Copyright (c) Tao Pang 1997.
C
      INTEGER I,N
      REAL H,G,C0,C1,C2,D,UTMP,Q(N),S(N),U(N)
C
      G = H*H/12.0
C
      DO      100  I = 2, N-1
        C0     = 1.0+G*Q(I-1)
        C1     = 2.0-10.0*G*Q(I)
        C2     = 1.0+G*Q(I+1)
        D      = G*(S(I+1)+S(I-1)+10.0*S(I))
        UTMP   = C1*U(I)-C0*U(I-1)+D
        U(I+1) = UTMP/C2
  100 CONTINUE
      RETURN
      END