ccccccccccccccccccccccccc     Program 3.3     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 SHOOTING
C
C Program for the boundary value problem with the shooting
C method.  The Runge-Kutta and secant methods are used.
C
      PARAMETER (N=101,M=5)
      DIMENSION Y(2,N)
C
C Initialization of the problem
C
      DL = 1.0E-06
      XL = 0.0
      XU = 1.0
      H  = (XU-XL)/(N-1)
      D  = 0.1
      YL = 0.0
      YU = 1.0
      X0 = (YU-YL)/(XU-XL)
      DX = 0.01
      X1 = X0+DX
C
C The secant search for the root
C
      Y(1,1) = YL
      DO 500 WHILE (ABS(D).GT.DL)
C
C The Runge-Kutta calculation of the first trial solution
C
        Y(2,1) = X0
        DO      100  I = 1, N-1
          X  = XL+H*I
          Y1 = Y(1,I)
          Y2 = Y(2,I)
          DK11 = H*G1(Y1,Y2,X)
          DK21 = H*G2(Y1,Y2,X)
          DK12 = H*G1((Y1+DK11/2.),(Y2+DK21/2.),(X+H/2.))
          DK22 = H*G2((Y1+DK11/2.),(Y2+DK21/2.),(X+H/2.))
          DK13 = H*G1((Y1+DK12/2.),(Y2+DK22/2.),(X+H/2.))
          DK23 = H*G2((Y1+DK12/2.),(Y2+DK22/2.),(X+H/2.))
          DK14 = H*G1((Y1+DK13),(Y2+DK23),(X+H))
          DK24 = H*G2((Y1+DK13),(Y2+DK23),(X+H))
          Y(1,I+1) = Y(1,I)+(DK11+2.*(DK12+DK13)+DK14)/6.
          Y(2,I+1) = Y(2,I)+(DK21+2.*(DK22+DK23)+DK24)/6.
  100   CONTINUE
        F0 = Y(1,N) - 1.0
C
C Runge-Kutta calculation of the second trial solution
C
        Y(2,1) = X1
        DO      200  I = 1, N-1
          X  = XL + H*I
          Y1 = Y(1,I)
          Y2 = Y(2,I)
          DK11 = H*G1(Y1,Y2,X)
          DK21 = H*G2(Y1,Y2,X)
          DK12 = H*G1((Y1+DK11/2.),(Y2+DK21/2.),(X+H/2.))
          DK22 = H*G2((Y1+DK11/2.),(Y2+DK21/2.),(X+H/2.))
          DK13 = H*G1((Y1+DK12/2.),(Y2+DK22/2.),(X+H/2.))
          DK23 = H*G2((Y1+DK12/2.),(Y2+DK22/2.),(X+H/2.))
          DK14 = H*G1((Y1+DK13),(Y2+DK23),(X+H))
          DK24 = H*G2((Y1+DK13),(Y2+DK23),(X+H))
          Y(1,I+1) = Y(1,I) + (DK11+2.*(DK12+DK13)+DK14)/6.
          Y(2,I+1) = Y(2,I) + (DK21+2.*(DK22+DK23)+DK24)/6.
  200   CONTINUE
        F1 = Y(1,N) - 1.0
C
        D = F1 - F0
        X2 = X1 - F1*(X1-X0)/D
        X0 = X1
        X1 = X2
  500 CONTINUE
      WRITE (6,999) (H*(I-1), Y(1,I),I=1,N,M)

C
  999 FORMAT (2F16.8)
      END
C
      FUNCTION G1(Y1,Y2,T)
        G1 = Y2
      RETURN
      END
C
      FUNCTION G2(Y1,Y2,T)
        PI = 4.0*ATAN(1.0)
        G2 =-PI*PI*(Y1+1.0)/4.0
      RETURN
      END