C***********************************************************************
C                                                                      *
C                   NONLINEAR SHOOTING METHOD                          *
C                                                                      *
C***********************************************************************
C
C
C
C     TO APPROXIMATE THE SOLUTION OF THE NONLINEAR BOUNDARY-VALUE
C     PROBLEM
C
C            Y'' = F(X,Y,Y'), A<=X<=B, Y(A)=ALPHA, Y(B)=BETA:
C
C
C     INPUT:   ENDPOINTS A,B; BOUNDARY CONDITIONS ALPHA, BETA; NUMBER OF
C              SUBINTERVALS N; TOLERANCE TOL; MAXIMUM NUMBER OF
C              ITERATIONS M.
C
C     OUTPUT:  APPROXIMATIONS W(1,I) TO Y(X(I)); W(2,I) TO Y'(X(I))
C              FOR EACH I=0,1,...,N OR A MESSAGE THAT THE MAXIMUM
C              NUMBER OF ITERATIONS WAS EXCEEDED.
C
C     INITIALIZATION
      DIMENSION W1(21),W2(21)
      CHARACTER NAME*14,NAME1*14,AA*1
      INTEGER INP,OUP,FLAG
      LOGICAL OK
      F(X,Y,Z)=(32+2*X**3-Y*Z)/8
C     FY(X,Y,Z) REPRESENTS PARTIAL OF F WITH RESPECT TO Y
      FY(X,Y,Z)=-Z/8
C     FYP(X,Y,Z) REPRESENTS PARTIAL OF F WITH RESPECT TO Y'
      FYP(X,Y,Z)=-Y/8
      OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      WRITE(6,*) 'This is the Nonlinear Shooting Method.'
      WRITE(6,*) 'Have the functions F, FY (partial of F'
      WRITE(6,*) 'with respect to y) and FYP (partial of'
      WRITE(6,*) 'F with respect to y'') been created? '
      WRITE(6,*) 'Enter Y or N '
      WRITE(6,*) ' '
      READ(5,*)  AA
      IF(( AA .EQ. 'Y' ) .OR. ( AA .EQ. 'y' )) THEN
         OK = .FALSE.
110      IF (OK) GOTO 11
            WRITE(6,*) 'Input left and right endpoints separated by'
            WRITE(6,*) 'blank - include decimal point'
            WRITE(6,*) ' '
            READ(5,*) A, B
            IF (A.GE.B) THEN
               WRITE(6,*) 'Left endpoint must be less'
               WRITE(6,*) 'than right endpoint'
            ELSE
               OK = .TRUE.
            ENDIF
         GOTO 110
11       OK = .FALSE.
         WRITE(6,*) 'Input Y(',A,')'
         WRITE(6,*) ' '
         READ(5,*) ALPHA
         WRITE(6,*) 'Input Y(',B,')'
         WRITE(6,*) ' '
         READ(5,*) BETA
12       IF (OK) GOTO 13
            WRITE(6,*) 'Input a positive integer for the number'
            WRITE(6,*) 'of subinvervals '
            WRITE(6,*) ' '
            READ(5,*) N
            IF ( N .LE. 0 ) THEN
              WRITE(6,*) 'Must be positive integer '
            ELSE
              OK = .TRUE.
            ENDIF
         GOTO 12
13       OK=.FALSE.
14       IF (OK) GOTO 15
            WRITE(6,*) 'Input tolerance '
            WRITE(6,*) ' '
            READ(5,*) TOL
            IF (TOL.LE.0.0) THEN
               WRITE(6,*) 'Tolerance must be positive '
               WRITE(6,*) ' '
            ELSE
               OK = .TRUE.
            ENDIF
         GOTO 14
15       OK = .FALSE.
16       IF (OK) GOTO 17
            WRITE(6,*) 'Input maximum number of iterations '
            WRITE(6,*) '- no decimal point '
            WRITE(6,*) ' '
            READ(5,*) M
            IF ( M .LE. 0 ) THEN
              WRITE(6,*) 'Must be positive integer '
              WRITE(6,*) ' '
            ELSE
              OK = .TRUE.
            ENDIF
         GOTO 16
17       CONTINUE
      ELSE
         WRITE(6,*) 'The program will end so that the functions'
         WRITE(6,*) 'F, FY and FYP can be created '
         OK = .FALSE.
      ENDIF
      IF(.NOT.OK) GOTO 400
      WRITE(6,*) 'Select output destinations: '
      WRITE(6,*) '1. Screen '
      WRITE(6,*) '2. Text file '
      WRITE(6,*) 'Enter 1 or 2 '
      WRITE(6,*) ' '
      READ(5,*) FLAG
      IF ( FLAG .EQ. 2 ) THEN
         WRITE(6,*) 'Input the file name in the form - '
         WRITE(6,*) 'drive:name.ext'
         WRITE(6,*) 'with the name contained within quotes'
         WRITE(6,*) 'as example:   ''A:OUTPUT.DTA'' '
         WRITE(6,*) ' '
         READ(5,*) NAME1
         OUP = 3
         OPEN(UNIT=OUP,FILE=NAME1,STATUS='NEW')
      ELSE
         OUP = 6
      ENDIF
      WRITE(OUP,*) 'NONLINEAR SHOOTING METHOD'
C     STEP 1
      H=(B-A)/N
      K=1
      TK=(BETA-ALPHA)/(B-A)
C     STEP 2
100   IF (K.GT.M) GOTO 200
C          STEP 3
C          W1(1)=W(1,0), W2(1)=W(2,0), U1=U(1), U2=U(2)
           W1(1)=ALPHA
           W2(1)= TK
           U1=0
           U2=1
C          STEP 4
C          RUNGE-KUTTA METHOD FOR SYSTEMS IS USED IN STEPS 5, 6
           DO 10 I=1,N
C               STEP 5
                X=A+(I-1)*H
                T=X+H/2
C               STEP 6
                XK11=H*W2(I)
                XK12=H*F(X,W1(I),W2(I))
                XK21=H*(W2(I)+XK12/2)
                XK22=H*F(T,W1(I)+XK11/2,W2(I)+XK12/2)
                XK31=H*(W2(I)+XK22/2)
                XK32=H*F(T,W1(I)+XK21/2,W2(I)+XK22/2)
                XK41=H*(W2(I)+XK32)
                XK42=H*F(X+H,W1(I)+XK31,W2(I)+XK32)
                W1(I+1)=W1(I)+(XK11+2*XK21+2*XK31+XK41)/6
                W2(I+1)=W2(I)+(XK12+2*XK22+2*XK32+XK42)/6
                YK11=H*U2
                YK12=H*(FY(X,W1(I),W2(I))*U1+FYP(X,W1(I),W2(I))*U2)
                YK21=H*(U2+YK12/2)
                YK22=H*(FY(T,W1(I),W2(I))*(U1+YK11/2)
     *          +FYP(T,W1(I),W2(I))*(U2+YK21/2))
                YK31=H*(U2+YK22/2)
                YK32=H*(FY(T,W1(I),W2(I))*(U1+YK21/2)
     *          +FYP(T,W1(I),W2(I))*(U2+YK22/2))
                YK41=H*(U2+YK32)
                YK42=H*(FY(X+H,W1(I),W2(I))*(U1+YK31)+FYP(X+H,W1(I),
     *          W2(I))*(U2+YK32))
                U1=U1+(YK11+2*YK21+2*YK31+YK41)/6
                U2=U2+(YK21+2*YK22+2*YK32+YK42)/6
10         CONTINUE
C          TEST FOR ACCURACY
C          STEP 7
           IF(ABS(W1(N+1)-BETA).LE.TOL) THEN
C               PROCESS IS COMPLETE
C               W1(J), W2(J) ARE THE APPROXIMATIONS FOR Y AND Y'
                WRITE(OUP,2) K
C               STEP 8
                DO 20 I=1,N
                     J=I+1
                     X=A+I*H
20              WRITE(OUP,1) X,W1(J),W2(J)
C               STEP 9
                GOTO 400
           END IF
C          STEP 10
C          NEWTON'S METHOD APPLIED TO IMPROVE TK
           TK=TK-(W1(N+1)-BETA)/U1
           K=K+1
      GOTO 100
C     STEP 11
C     METHOD FAILED
200   WRITE(OUP,3) M
400   CLOSE(UNIT=5)
      CLOSE(UNIT=OUP)
      IF(OUP.NE.6) CLOSE(UNIT=6)
      STOP
1     FORMAT(1X,3(E15.8,3X))
2     FORMAT(1X,'ORDER OF OUTPUT - X(I),W1(I),W2(I)',1X,I3,1X,'ITERATION
     *S')
3     FORMAT(1X,'METHOD FAILED AFTER ITERATION NO.' ,I4)
      END
