ccccccccccccccccccccccccc     Program 3.C     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 WAVE
C
C Program for the eigenvalue problem with a combination of the
C bisection method and the Numerov algorithm for u" = -k**2*u
C with boundary conditions u(0)=u(1)=0.
C Copyright (c) Tao Pang 1997.
C
      PARAMETER (N=101)
      INTEGER ISTEP
      REAL H,AK,BK,DK,EK,DL,F0,F1,Q(N),S(N),U(N)
C
C Initialization of the problem
C
      H     = 1.0/(N-1)
      AK    = 2.0
      BK    = 4.0
      DK    = 1.0
      DL    = 1.0E-06
      ISTEP = 0
      U(1)  = 0.0
      U(2)  = 0.01
      EK    = AK
      DO    100  I = 1,N
        S(I) = 0.0
        Q(I) = EK*EK
  100 CONTINUE
      CALL NMRV (N,H,Q,S,U)
      F0 = U(N)
C
C Bisection method for the root
C
      DO      200  WHILE (ABS(DK).GT.DL)
        EK = (AK+BK)/2.0
        DO    190  I = 1,N
          Q(I) = EK*EK
  190   CONTINUE
        CALL NMRV (N,H,Q,S,U)
        F1 = U(N)
        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
      WRITE (6,999) ISTEP,EK,DK,F1
  999 FORMAT (I4,3F16.8)
      STOP
      END
C
      SUBROUTINE NMRV (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.77)-(3.80) 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)-Q(I+1))/2.0+Q(I))
        C1     = 2.0-G*(Q(I+1)+Q(I-1)+8.0*Q(I))
        C2     = 1.0+G*((Q(I+1)-Q(I-1))/2.0+Q(I))
        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