C***********************************************************************
C                                                                      *
C     ADAMS VARIABLE STEP-SIZE PREDICTOR-CORRECTOR METHOD              *
C                                                                      *
C***********************************************************************
C
C
C
C     TO APPROXIMATE THE SOLUTION OF THE INITIAL VALUE PROBLEM
C                Y'=F(T,Y), A <= T <= B , Y(A)=ALPHA,
C     WITH LOCAL TRUNCATION ERROR WITHIN A GIVEN TOLERANCE:
C
C     INPUT:   ENDPOINTS A,B; INITIAL CONDITION ALPHA; TOLERANCE
C              TOL; MAXIMUM STEPSIZE HMAX; MINIMUM STEPSIZE HMIN.
C
C     OUTPUT:  I,T(I),W(I),H WHERE AT THE ITH STEP W(I) APPROXIMATES
C              Y(T(I)) AND STEPSIZE H WAS USED OR A MESSAGE THAT THE
C              MINIMUM STEP SIZE WAS EXCEEDED.
C
C     STEP 1 SETS UP THE SUBALGORITHM FOR RK4
C     THE FUNCTION F MUST ALSO BE DEFINED IN THE SUBROUTINE
C
      DIMENSION T(100),W(100)
      CHARACTER NAME*14,NAME1*14,AA*1
      INTEGER INP,OUP,FLAG
      LOGICAL OK
      F(TZ,WZ)= WZ-TZ*TZ+1
      OPEN(UNIT=5,FILE='CON',ACCESS='SEQUENTIAL')
      OPEN(UNIT=6,FILE='CON',ACCESS='SEQUENTIAL')
      WRITE(6,*) 'This is the Adams Variable Step-size'
      WRITE(6,*) 'Predictor-Corrector Method.'
      WRITE(6,*) 'Has the function F been created in the program? '
      WRITE(6,*) 'NOTE:  The function is defined in two places.'
      WRITE(6,*) 'Enter Y or N '
      WRITE(6,*) ' '
      READ(5,*)  AA
      IF(( AA .EQ. 'Y' ) .OR. ( AA .EQ. 'y' )) THEN
         OK = .FALSE.
10       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 10
11       OK = .FALSE.
13       WRITE(6,*) 'Input the initial condition.'
         WRITE(6,*) ' '
         READ(5,*) ALPHA
         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 '
            ELSE
               OK = .TRUE.
            ENDIF
         GOTO 14
15       OK = .FALSE.
16       IF (OK) GOTO 17
            WRITE(6,*) 'Input minimum and maximum mesh spacing'
            WRITE(6,*) 'separated by blank.'
            WRITE(6,*) ' '
            READ(5,*) HMIN, HMAX
            IF ((HMIN.GE.HMAX) .OR. (HMIN.LE.0.0)) THEN
               WRITE(6,*) 'Minimum mesh spacing must be positive '
               WRITE(6,*) 'and less than maximum mesh spacing.'
            ELSE
               OK = .TRUE.
            ENDIF
         GOTO 16
17       CONTINUE
      ELSE
         WRITE(6,*) 'The program will end so that the function F '
         WRITE(6,*) '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,*) 'ADAMS VARIABLE PREDICTOR-CORRECTOR METHOD'
      WRITE(OUP,*) 'The fifth column is sigma.  Output only if'
      WRITE(OUP,*) 'approximation is from Adams formulas.'
      WRITE(OUP,6)
6     FORMAT(3X,'I',12X,'t(i)',12X,'w(i)',15x,'h')
C     STEP 2
C     SUBSCRIPTS ARE SHIFTED TO AVOID ZERO SUBSCRIPTS
      T(1)= A
      W(1)= ALPHA
      H= HMAX
C     FLAG WILL BE USED TO EXIT THE LOOP IN STEP 4
      FLAG=1
C     LAST WILL INDICATE WHEN THE LAST VALUE IS CALCULATED.
      LAST=0
      WRITE(OUP,1) T(1),W(1)
C     STEP 3
      KK=1
      CALL RK4(KK,H,W,T)
C     NFLAG INDICATES COMPUTATION FROM RK4
      NFLAG=1
      I=5
C     USE TT IN PLACE OF T
      TT=T(4)+H
C     STEP 4
100   IF (FLAG .NE. 1 ) GOTO 200
C         STEP 5
C         PREDICT W(I)
          WP=W(I-1)+H*(55*F(T(I-1),W(I-1))-59*F(T(I-2),W(I-2))+37*F(T(I-
     1    3),W(I-3))-9*F(T(I-4),W(I-4)))/24
C         CORRECT W(I)
          WC=W(I-1)+H*(9*F(TT,WP)+19*F(T(I-1),W(I-1))-5*F(T(I-2),W(I-2))
     1    +F(T(I-3),W(I-3)))/24
          SIG=19*(ABS(WC-WP))/(270*H)
C         STEP 6
          IF(SIG.LE.TOL) THEN
C             STEP 7
C             RESULT ACCEPTED
              W(I)=WC
              T(I)=TT
C             STEP 8
              IF(NFLAG.EQ.1) THEN
                   K=I-3
                   KK=I-1
                   DO 20 J=K,KK
                   JJ=J-1
20                 WRITE(OUP,2) JJ,T(J),W(J),H
                   J=I-1
                   WRITE(OUP,3) J,T(I),W(I),H,SIG
              ELSE
                   J=I-1
                   WRITE(OUP,3) J,T(I),W(I),H,SIG
              END IF
C             STEP 9
              IF (LAST.EQ.1) THEN
                  FLAG=0
              ELSE
C                  STEP 10
                   I=I+1
                   NFLAG=0
C                  STEP 11
                   IF((SIG.LE.0.1*TOL).OR.(T(I-1)+H.GT.B)) THEN
C                  INCREASE H IF IT MORE ACCURATE THAN REQUIRED OR
C                  DECREASE H TO INCLUDE B AS A MESH POINT.
C                  STEP 12
C                  TO AVOID UNDERFLOW
                        IF (SIG .LE. 1.0E-70) THEN
                             Q=4
                        ELSE
                             Q=(TOL/(2*SIG))**.25
                        END IF
C                       STEP 13
                        IF(Q.GT.4) THEN
                             H=4*H
                        ELSE
                             H=Q*H
                        END IF
C                       STEP 14
                        IF(H.GT.HMAX) H=HMAX
C                       STEP 15
C                  AVOID TERMINATING WITH CHANGE IN STEPSIZE
                        IF(T(I-1)+4*H.GE.B) THEN
                            H=(B-T(I-1))/4
                            LAST=1
                        ENDIF
C                       STEP 16
                        NFLAG=1
                        KK=I-1
                        CALL RK4(KK,H,W,T)
                        I=I+3
                   END IF
              END IF
          ELSE
C             RESULT REJECTED
C             STEP 17
              Q=(TOL/(2*SIG))**.25
C             STEP 18
              IF(Q.LT.0.1) THEN
                   H=0.1*H
              ELSE
                   H=Q*H
              END IF
C             STEP 19
              IF(H.LT.HMIN) THEN
C                  PROCEDURE FAILS
                   WRITE(OUP,4)
                   FLAG=0
              ELSE
C             PREVIOUS RESULTS ALSO REJECTED
                   IF(NFLAG.EQ.1) I=I-3
                   KK=I-1
                   CALL RK4(KK,H,W,T)
                   I=I+3
                   NFLAG=1
              ENDIF
          END IF
C     STEP 20
      TT=T(I-1)+H
      GOTO 100
C     STEP 21
200   CONTINUE
400   CLOSE(UNIT=5)
      CLOSE(UNIT=OUP)
      IF(OUP.NE.6) CLOSE(UNIT=6)
      STOP
1     FORMAT(5X,E15.8,1X,E15.8)
2     FORMAT(1X,I3,1X,E15.8,1X,E15.8,1X,E15.8)
3     FORMAT(1X,I3,1X,E15.8,1X,E15.8,1X,E15.8,1X,E15.8)
4     FORMAT(1X,'HMIN EXCEEDED - PROCEDURE FAILS')
      END
      SUBROUTINE RK4(K,H,V,X)
      DIMENSION V(100),X(100)
      F(XZ,VZ)=VZ-XZ*XZ+1
      DO 20 I=1,3
      J=K+I-1
      XK1=H*F(X(J),V(J))
      XK2=H*F(X(J)+H/2,V(J)+XK1/2)
      XK3=H*F(X(J)+H/2,V(J)+XK2/2)
      XK4=H*F(X(J)+H,V(J)+XK3)
      V(J+1)=V(J)+(XK1+2*(XK2+XK3)+XK4)/6
20    X(J+1)=X(J)+H
      RETURN
      END
