      PROGRAM SOLVEDE
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 9.1 (Euler's Method).
C     Section 9.2, Euler's Method, Page 435
C
C     Algorithm 9.2 (Heun's Method).
C     Section 9.3, Heun's Method, Page 441
C
C     Algorithm 9.3 (Taylor's Method of Order 4).
C     Section 9.4, Taylor Series Method, Page 448

C     Algorithm 9.4 (Runge-Kutta Method of Order 4).
C     Section 9.5, Runge-Kutta Methods, Page 460
C
      PARAMETER(MaxM=1000)
      INTEGER DoMo,M,Mend,Meth,Order,State
      REAL A,B,D,T,TK,Y,YK,Y0
      DIMENSION D(1:4),T(0:1000),Y(0:1000)
      CHARACTER Ans*1
      EXTERNAL DF,F
      Meth=1
      A=0
      B=0
      Y0=0
      M=1
      State=1
10    CONTINUE
      CALL MESSAGE(Meth,Order,State)
      DoMo=1
20      CONTINUE
          CALL INPUT(Meth)
30        CONTINUE
            CALL EPOINTS(A,B,Y0,M,State)
            IF (Meth.EQ.1) THEN
              CALL EULER(F,A,B,Y0,M,Mend,T,Y)
            ELSEIF (Meth.EQ.2) THEN
              CALL CAUCHY(F,A,B,Y0,M,Mend,T,Y)
            ELSEIF (Meth.EQ.3) THEN
              CALL HEUN(F,A,B,Y0,M,Mend,T,Y)
            ELSEIF (Meth.EQ.4) THEN
              CALL TAYLOR(DF,A,B,Y0,M,Mend,T,TK,Y,YK,D,Order)
            ELSEIF (Meth.EQ.5) THEN
              CALL TAYLOR(DF,A,B,Y0,M,Mend,T,TK,Y,YK,D,Order)
            ELSEIF (Meth.EQ.6) THEN
              CALL RUNGE(F,A,B,Y0,M,Mend,T,Y)
            ENDIF
            CALL RESULTS(T,Y,M,Mend,Meth)
            WRITE(9,*)' '
            WRITE(9,*)'      Want to try another initial condition ?
     & <Y/N> '
            READ(9,'(A)') Ans
            IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN
              State=2
            ELSE
              State=0
            ENDIF
          IF (State.EQ.0.OR.State.EQ.1) GOTO 30
          WRITE(9,*)'   Want to change the differential equation ?
     & <Y/N> '
            READ(9,'(A)') Ans
          IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN
            DoMo=0
          ELSE
            State=0
          ENDIF
        IF (DoMo.EQ.1) GOTO 20
      WRITE(9,*)'Want to try another method of approximation ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN
        Meth=0
      ELSE
        State=0
      ENDIF
      IF (Meth.NE.0) GOTO 10
      RETURN
      END

      REAL FUNCTION F(T,Y)
      REAL T,Y
      F=T**2+Y**2
      RETURN
      END

      SUBROUTINE DF(T,Y,D)
      REAL D,T,Y
      DIMENSION D(1:4)
        D(1) = ( T - Y )/2
        D(2) = ( 2 - T + Y)/4
        D(3) = (-2 + T - Y)/8
        D(4) = ( 2 - T + Y)/16
        RETURN
      END

      SUBROUTINE PRINTFUN
        WRITE(9,*)'          Y` = T**2+Y**2'
      RETURN
      END

      SUBROUTINE EULER(F,A,B,Y0,M,Mend,T,Y)
      PARAMETER(Big=1E13)
      INTEGER K,M,Mend
      REAL A,B,H,T,Y,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO K=0,M-1
        Y(K+1)=Y(K)+H*F(T(K),Y(K))
        T(K+1)=A+H*(K+1)
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) RETURN
      ENDDO
      RETURN
      END

      SUBROUTINE HEUN(F,A,B,Y0,M,Mend,T,Y)
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend
      REAL A,B,H,K1,K2,P,T,Y,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO J=0,M-1
        K1=F(T(J),Y(J))
        P=Y(J)+H*K1
        T(J+1)=A+H*(J+1)
        Mend=J+1
        IF (Big.LT.ABS(P)) THEN
          Y(J+1)=P
          RETURN
        ENDIF
        K2=F(T(J+1),P)
        Y(J+1)=Y(J)+H*(K1+K2)/2
        IF (Big.LT.ABS(Y(J+1))) RETURN
      ENDDO
      RETURN
      END

      SUBROUTINE CAUCHY(F,A,B,Y0,M,Mend,T,Y)
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend
      REAL A,B,H,K1,K2,P,T,TP,Y,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO J=0,M-1
        K1=F(T(J),Y(J))
        P=Y(J)+H*K1/2
        T(J+1)=A+H*(J+1)
        Mend=J+1
        IF (Big.LT.ABS(P)) THEN
          Y(J+1)=P
          RETURN
        ENDIF
        TP=T(J)+H/2
        K2=F(TP,P)
        Y(J+1)=Y(J)+H*K2
        IF (Big.LT.ABS(Y(J+1))) RETURN
      ENDDO
      RETURN
      END

      SUBROUTINE TAYLOR(DF,A,B,Y0,M,Mend,T,TK,Y,YK,D,Order)
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend,Order
      REAL A,B,D,H,T,TK,Y,YK,Y0
      DIMENSION D(1:4),T(0:1000),Y(0:1000)
      EXTERNAL DF
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO K=0,M-1
        TK=T(K)
        YK=Y(K)
        CALL DF(TK,YK,D)
        IF (Order.EQ.3) THEN
          Y(K+1)=YK+H*(D(1)+H*(D(2)/2+H*D(3)/6))
        ELSE
          Y(K+1)=YK+H*(D(1)+H*(D(2)/2+H*(D(3)/6+H*D(4)/24)))
        ENDIF
        T(K+1)=A+H*(K+1)
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) RETURN
      ENDDO
      RETURN
      END

      SUBROUTINE RUNGE(F,A,B,Y0,M,Mend,T,Y)
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend,Pole
      REAL A,B,H,K1,K2,K3,K4,P,T,TJ,Y,YJ,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO J=0,M-1
        TJ=T(J)
        YJ=Y(J)
        K1=H*F(TJ,YJ)
          IF (Big.LT.ABS(YJ+0.5*K1)) RETURN
        K2=H*F(TJ+H/2,YJ+0.5*K1)
          IF (Big.LT.ABS(YJ+0.5*K2)) RETURN
        K3=H*F(TJ+H/2,YJ+0.5*K2)
          IF (Big.LT.ABS(YJ+K3)) RETURN
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
        Mend=J+1
          IF (Big.LT.ABS(Y(J+1))) RETURN
      ENDDO
      RETURN
      END

      SUBROUTINE XRUNGE(F,A,B,Y0,M,Mend,T,Y)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend,Pole
      REAL A,B,H,K1,K2,K3,K4,P,T,TJ,Y,YJ,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      Pole=0
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO J=0,M-1
        TJ=T(J)
        YJ=Y(J)
        Mend=J+1
        K1=H*F(TJ,YJ)
        IF (Big.LT.ABS(YJ+0.5*K1)) THEN
          Pole=1
          GOTO 10
        ENDIF
        K2=H*F(TJ+H/2,YJ+0.5*K1)
        IF (Big.LT.ABS(YJ+0.5*K2)) THEN
          Pole=1
          GOTO 10
        ENDIF
        K3=H*F(TJ+H/2,YJ+0.5*K2)
        IF (Big.LT.ABS(YJ+K3)) THEN
          Pole=1
          GOTO 10
        ENDIF
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
        IF (Big.LT.ABS(Y(J+1))) GOTO 10
      ENDDO
10    CONTINUE
      IF (Pole.EQ.0) RETURN
      P=YJ+K1
      T(J+1)=A+H*(J+1)
      IF (Big.LT.ABS(P)) THEN
        Y(J+1)=P
        RETURN
      ENDIF
      K2=F(T(J+1),P)
      Y(J+1)=YJ+H*(K1+K2)/2
      RETURN
      END

      SUBROUTINE XEULER(F,A,B,Y0,M,Mend,T,Y)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      INTEGER K,M,Mend
      REAL A,B,H,T,Y,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO 10 K=0,M-1
        Y(K+1)=Y(K)+H*F(T(K),Y(K))
        T(K+1)=A+H*(K+1)
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) GOTO 20
10    CONTINUE
20    CONTINUE
      RETURN
      END

      SUBROUTINE XHEUN(F,A,B,Y0,M,Mend,T,Y)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend
      REAL A,B,H,K1,K2,P,T,Y,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO 10 J=0,M-1
      K1=F(T(J),Y(J))
      P=Y(J)+H*K1
      T(J+1)=A+H*(J+1)
      Mend=J+1
      IF (Big.LT.ABS(P)) THEN
        Y(J+1)=P
        GOTO 20
      ENDIF
      K2=F(T(J+1),P)
      Y(J+1)=Y(J)+H*(K1+K2)/2
      IF (Big.LT.ABS(Y(J+1))) GOTO 20
10    CONTINUE
20    CONTINUE
      RETURN
      END

      SUBROUTINE XCAUCHY(F,A,B,Y0,M,Mend,T,Y)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend
      REAL A,B,H,K1,K2,P,T,TP,Y,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO 10 J=0,M-1
      K1=F(T(J),Y(J))
      P=Y(J)+H*K1/2
      T(J+1)=A+H*(J+1)
      Mend=J+1
      IF (Big.LT.ABS(P)) THEN
        Y(J+1)=P
        GOTO 20
      ENDIF
      TP=T(J)+H/2
      K2=F(TP,P)
      Y(J+1)=Y(J)+H*K2
      IF (Big.LT.ABS(Y(J+1))) GOTO 20
10    CONTINUE
20    CONTINUE
      RETURN
      END

      SUBROUTINE XTAYLOR(DF,A,B,Y0,M,Mend,T,TK,Y,YK,D,Order)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend,Order
      REAL A,B,D,H,T,TK,Y,YK,Y0
      DIMENSION D(1:4),T(0:1000),Y(0:1000)
      EXTERNAL DF
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO 10 K=0,M-1
      TK=T(K)
      YK=Y(K)
      CALL DF(TK,YK,D)
      IF (Order.EQ.3) THEN
        Y(K+1)=YK+H*(D(1)+H*(D(2)/2+H*D(3)/6))
      ELSE
        Y(K+1)=YK+H*(D(1)+H*(D(2)/2+H*(D(3)/6+H*D(4)/24)))
      ENDIF
      T(K+1)=A+H*(K+1)
      Mend=K+1
      IF (Big.LT.ABS(Y(K+1))) GOTO 20
10    CONTINUE
20    CONTINUE
      RETURN
      END

      SUBROUTINE XRUNGE(F,A,B,Y0,M,Mend,T,Y)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      INTEGER J,M,Mend,Pole
      REAL A,B,H,K1,K2,K3,K4,P,T,TJ,Y,YJ,Y0
      DIMENSION T(0:1000),Y(0:1000)
      EXTERNAL F
      Pole=0
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      DO 10 J=0,M-1
        TJ=T(J)
        YJ=Y(J)
        Mend=J+1
        K1=H*F(TJ,YJ)
        IF (Big.LT.ABS(YJ+0.5*K1)) THEN
          Pole=1
          GOTO 20
        ENDIF
        K2=H*F(TJ+H/2,YJ+0.5*K1)
        IF (Big.LT.ABS(YJ+0.5*K2)) THEN
          Pole=1
          GOTO 20
        ENDIF
        K3=H*F(TJ+H/2,YJ+0.5*K2)
        IF (Big.LT.ABS(YJ+K3)) THEN
          Pole=1
          GOTO 20
        ENDIF
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
        IF (Big.LT.ABS(Y(J+1))) GOTO 20
10    CONTINUE
20    CONTINUE
      IF (Pole.EQ.0) GOTO 30
      P=YJ+K1
      T(J+1)=A+H*(J+1)
      IF (Big.LT.ABS(P)) THEN
        Y(J+1)=P
        GOTO 30
      ENDIF
      K2=F(T(J+1),P)
      Y(J+1)=YJ+H*(K1+K2)/2
30    CONTINUE
      RETURN
      END

      SUBROUTINE MESSAGE(Meth,Order,State)
      INTEGER K,Meth,Order,State
      REAL Resp
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'          Solution of the differential equation  Y` = F
     +(T,Y)'
      WRITE(9,*)' '
      WRITE(9,*)'          with the initial condition  Y(A) = Y .'
      WRITE(9,*)'                                              0'
      WRITE(9,*)' '
      WRITE(9,*)'          A numerical approximation is computed over [A
     +,B].'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'          Choose the method of approximation:'
      WRITE(9,*)' '
      WRITE(9,*)'    < 1 > Euler`s method'
      WRITE(9,*)' '
      WRITE(9,*)'    < 2 > Modified Euler-Cauchy method'
      WRITE(9,*)' '
      WRITE(9,*)'    < 3 > Heun`s method'
      WRITE(9,*)' '
      WRITE(9,*)'    < 4 > Taylor`s method of order N=3'
      WRITE(9,*)' '
      WRITE(9,*)'    < 5 > Taylor`s method of order N=4'
      WRITE(9,*)' '
      WRITE(9,*)'    < 6 > Runge-Kutta method of order N=4'
      WRITE(9,*)' '
      WRITE(9,*)'          SELECT < 1 - 6 > ? '
      Resp=Meth
      READ(9,*) Resp
      Meth=INT(Resp)
      IF ((Meth.LT.1.OR.Meth.GT.6).AND.State.NE.0) THEN
        Meth=1
      ENDIF
      IF (Meth.EQ.4) THEN
        Order=3
      ENDIF
      IF (Meth.EQ.5) THEN
        Order=4
      ENDIF
      RETURN
      END

      SUBROUTINE INPUT(Meth)
      INTEGER Meth
      REAL Valu
      CHARACTER Ans*1
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'     You chose '
      WRITE(9,*)' '
      IF (Meth.EQ.1) THEN
        WRITE(9,*)'     Euler`s method to solve the D.E.'
      ELSEIF (Meth.EQ.2) THEN
        WRITE(9,*)'     the Modified Euler-Cauchy method to solve
     + the D.E.'
      ELSEIF (Meth.EQ.3) THEN
        WRITE(9,*)'     Heun`s method to solve the D.E.'
      ELSEIF (Meth.EQ.4) THEN
        WRITE(9,*)'     Taylor`s method of order N=3 to solve
     + the D.E.'
      ELSEIF (Meth.EQ.5) THEN
        WRITE(9,*)'     Taylor`s method of order N=4 to solve
     + the D.E.'
      ELSEIF (Meth.EQ.6) THEN
        WRITE(9,*)'     the Runge-Kutta method of order N=4 to solve
     + the D.E.'
      ENDIF
      WRITE(9,*)' '
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)'     with the initial condition Y(A) = Y0.'
      WRITE(9,*)' '
      WRITE(9,*)'     A numerical approximation is computed over [A,B].'
      WRITE(9,*)' '
      WRITE(9,*)'     You must enter the endpoints for the interval,'
      WRITE(9,*)' '
      WRITE(9,*)'     the initial condition Y0, and the number of steps
     + M.'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'         Press the <ENTER> key. '
      READ(9,'(A)') Ans
      RETURN
      END

      SUBROUTINE EPOINTS(A,B,Y0,M,State)
      INTEGER I,M,Stat,State
      REAL A,B,Y0,Valu
C     STATUS = (Change,Enter,Done)
      CHARACTER Ans*1
      Stat=1
      IF (State.EQ.0) THEN
        Stat=0
      ENDIF
10    CONTINUE
      DO 20 I=1,15
        WRITE(9,*)' '
20    CONTINUE
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)' '
      IF (Stat.EQ.1) THEN
        WRITE(9,*)'     ENTER  the  left  endpoint  A = '
        READ(9,*) A
        WRITE(9,*)' '
        WRITE(9,*)'     ENTER  the  right endpoint  B = '
        READ(9,*) B
        WRITE(9,*)' '
        WRITE(9,*)'     ENTER initial condition  Y(A) = '
        READ(9,*) Y0
        WRITE(9,*)' '
        WRITE(9,*)'     ENTER the number of steps   M = '
        Valu=M
        READ(9,*) Valu
        M=INT(Valu)
        IF (M.LT.1) THEN
          M=1
        ENDIF
        IF (M.GT.1000) THEN
          M=1000
        ENDIF
      ELSE
        WRITE(9,*)'     The left  endpoint is    A =',A
        WRITE(9,*)' '
        WRITE(9,*)'     The right endpoint is    B =',B
        WRITE(9,*)' '
        WRITE(9,*)'     Initial condition is  Y(A) =',Y0
        WRITE(9,*)' '
        WRITE(9,*)'     The number of steps is   M = ',M
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'     Do you want to make a change ?  <Y/N>  '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y'.OR.Ans.EQ.'y') THEN
          Stat=0
          DO 30 I=1,15
            WRITE(9,*)' '
30        CONTINUE
          WRITE(9,*)'          '
          CALL PRINTFUN
          WRITE(9,*)' '
          WRITE(9,*)' '
          WRITE(9,*)'     The current left  endpoint is A =',A
            WRITE(9,*)'     ENTER  the NEW left  endpoint A = '
          READ(9,*) A
          WRITE(9,*)' '
          WRITE(9,*)'     The current right endpoint is B =',B
            WRITE(9,*)'     ENTER  the NEW right endpoint B = '
          READ(9,*) B
          WRITE(9,*)' '
          WRITE(9,*)'        The current   I. C. is  Y(A) =',Y0
            WRITE(9,*)'        ENTER the NEW    I. C.  Y(A) = '
          READ(9,*) Y0
          WRITE(9,*)' '
          WRITE(9,*)'        The current value of M is  M = ',M
            WRITE(9,*)'        ENTER  the  NEW  value of  M = '
          Valu=M
          READ(9,*) Valu
          M=INT(Valu)
          IF (M.LT.1) THEN
            M=1
          ENDIF
          IF (M.GT.1000) THEN
            M=1000
          ENDIF
      ELSE
        Stat=2
      ENDIF
      IF (Stat.EQ.0.OR.Stat.EQ.1) GOTO 10
      RETURN
      END

      SUBROUTINE RESULTS(T,Y,M,Mend,Meth)
      INTEGER I,K,M,Mend,Meth
      REAL T,Y
      DIMENSION T(0:1000),Y(0:1000)
      DO 10 I=1,10
        WRITE(9,*)' '
10    CONTINUE
      IF (Meth.EQ.1) THEN
        WRITE(9,*)'Euler`s method was used to solve the D.E.'
      ELSEIF (Meth.EQ.2) THEN
        WRITE(9,*)'The Modified Euler-Cauchy method was used to solve th
     +e D.E.'
      ELSEIF (Meth.EQ.3) THEN
        WRITE(9,*)'Heun`s method was used to solve the D.E.'
      ELSEIF (Meth.EQ.4) THEN
        WRITE(9,*)'Taylor`s method of order N=3 was used to solve the D.
     +E.'
      ELSEIF (Meth.EQ.5) THEN
        WRITE(9,*)'Taylor`s method of order N=4 was used to solve the D.
     +E.'
      ELSEIF (Meth.EQ.6) THEN
        WRITE(9,*)'The Runge-Kutta method of order N=4 was used to solve
     + the D.E.'
      ENDIF
      WRITE(9,*)' '
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)'with  Y(',T(0),'  ) =',Y(0)
      WRITE(9,*)' '
      WRITE(9,*)'   K            T(K)           Y(K)'
      WRITE(9,*)' ---------------------------------------'
      DO 20 K=0,Mend
        WRITE(9,999) K, T(K), Y(K)
20    CONTINUE
      IF (Mend.LT.M) THEN
        WRITE(9,*)' '
        WRITE(9,*)'The solution points are approaching a pole.'
      ENDIF
999   FORMAT(I4,4X,4F15.7)
      PAUSE
      RETURN
      END

      PROGRAM SOLVEDE
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 9.5 (Runge-Kutta-Fehlberg Method RKF45).
C     Section 9.5, Runge-Kutta Methods, Page 461
C
C     Algorithm 9.6 (Adams-Bashforth-Moulton Method).
C     Section 9.6, Predictor-Corrector Method, Page 471
C
C     Algorithm 9.7 (Milne-Simpson Method).
C     Section 9.6, Predictor-Corrector Method, Page 472

C     Algorithm 9.8 (The Hamming Method).
C     Section 9.6, Predictor-Corrector Method, Page 473
C
      PARAMETER(MaxM=1000)
      INTEGER DoMo,M,Mend,Meth,Order,State
      REAL A,B,D,T,TK,Tol,Y,YK,Y0
      DIMENSION D(1:4),T(1:1000),Y(1:1000)
      CHARACTER Ans*1
      EXTERNAL DF,F
      Meth=1
      A=0
      B=0
      Y0=0
      M=1
      State=1
10    CONTINUE
      CALL MESSAGE(Meth,Order,State)
      DoMo=1
20      CONTINUE
          CALL INPUT(Meth)
30        CONTINUE
            CALL EPOINTS(A,B,Y0,M,Tol,State)
            IF (Meth.EQ.1) THEN
              CALL ABM(F,A,B,Y0,M,T,Y,Mend)
            ELSEIF (Meth.EQ.2) THEN
              CALL MILNE(F,A,B,Y0,M,T,Y,Mend)
            ELSEIF (Meth.EQ.3) THEN
              CALL HAMMING(F,A,B,Y0,M,T,Y,Mend)
            ELSEIF (Meth.EQ.6) THEN
              CALL RKF45(F,A,B,Y0,Tol,T,Y,M,MaxM,Mend)
            ENDIF
            CALL RESULTS(T,Y,M,Mend,Meth)
            WRITE(9,*)' '
            WRITE(9,*)'      Want to try another initial condition ? <Y/
     +N> '
            READ(9,'(A)') Ans
            IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN
              State=2
            ELSE
              State=0
            ENDIF
          IF (State.EQ.0.OR.State.EQ.1) GOTO 30
          WRITE(9,*)'   Want to change the differential equation ? <Y/N>
     + '
            READ(9,'(A)') Ans
          IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN
            DoMo=0
          ELSE
            State=0
          ENDIF
        IF (DoMo.EQ.1) GOTO 20
      WRITE(9,*)'Want to try another method of approximation ? <Y/N> '
      READ(9,'(A)') Ans
      IF (Ans.NE.'Y'.AND.Ans.NE.'y') THEN
        Meth=0
      ELSE
        State=0
      ENDIF
      IF (Meth.NE.0) GOTO 10
      RETURN
      END

      REAL FUNCTION F(T,Y)
      REAL T,Y
        F=T**2+Y**2
        RETURN
      END

      SUBROUTINE PRINTFUN
        WRITE(9,*)'          Y` = T**2+Y**2'
        RETURN
      END

      SUBROUTINE HAMMING(F,A,B,Y0,M,T,Y,Mend)
      PARAMETER(Big=1E13)
      REAL A,B,Cnew,Cold,H,F1,F2,F3,F4,Pmod,Pnew,Pold,T,TK,Y,YK,Y0
      REAL K1,K2,K3,K4,TJ,YJ
      DIMENSION T(1:1000),Y(1:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      !Use Runge-Kutta to compute three starting values.
      DO J=0,2
        TJ=T(J)
        YJ=Y(J)
        K1=H*F(TJ,YJ)
        K2=H*F(TJ+H/2,YJ+0.5*K1)
        K3=H*F(TJ+H/2,YJ+0.5*K2)
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
      ENDDO
      !Start of Procedure Hamming
      F1=F(T(1),Y(1))
      F2=F(T(2),Y(2))
      F3=F(T(3),Y(3))
      Pold=0
      Cold=0
      DO K=3,M-1
        Pnew=Y(K-3)+4*H*(2*F1-F2+2*F3)/3
        Pmod=Pnew+112*(Cold-Pold)/121
        T(K+1)=A+H*(K+1)
        F4=F(T(K+1),Pmod)
        Cnew=(9*Y(K)-Y(K-2)+3*H*(-F2+2*F3+F4))/8
        Y(K+1)=Cnew+9*(Pnew-Cnew)/121
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) RETURN
        Pold=Pnew
        Cold=Cnew
        F1=F2
        F2=F3
        F3=F(T(K+1),Y(K+1))
      ENDDO
      RETURN
      END

      SUBROUTINE ABM(F,A,B,Y0,M,T,Y,Mend)
      PARAMETER(Big=1E13)
      INTEGER J,K,M
      REAL A,B,H,H2,F0,F1,F2,F3,F4,P,T,TK,Y,YK,Y0
      REAL K1,K2,K3,K4,TJ,YJ
      DIMENSION T(1:1000),Y(1:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      !Use Runge-Kutta to compute three starting values.
      DO J=0,2
        TJ=T(J)
        YJ=Y(J)
        K1=H*F(TJ,YJ)
        K2=H*F(TJ+H/2,YJ+0.5*K1)
        K3=H*F(TJ+H/2,YJ+0.5*K2)
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
      ENDDO
      !Start of Procedure Adams-Bashforth-Moulton
      F0=F(T(0),Y(0))
      F1=F(T(1),Y(1))
      F2=F(T(2),Y(2))
      F3=F(T(3),Y(3))
      H2=H/24
      DO K=3,M-1
        TK=T(K)
        YK=Y(K)
        P=YK+H2*(-9*F0+37*F1-59*F2+55*F3)
        T(K+1)=A+H*(K+1)
        F4=F(T(K+1),P)
        Y(K+1)=YK+H2*(F1-5*F2+19*F3+9*F4)
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) RETURN
        F0=F1
        F1=F2
        F2=F3
        F3=F(T(K+1),Y(K+1))
      ENDDO
      RETURN
      END

      SUBROUTINE MILNE(F,A,B,Y0,M,T,Y,Mend)
      PARAMETER(Big=1E13)
      INTEGER J,K,M
      REAL A,B,H,F1,F2,F3,F4,Pmod,Pnew,Pold,T,TK,Y,YK,Yold,Y0
      REAL K1,K2,K3,K4,TJ,YJ
      DIMENSION T(1:1000),Y(1:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
      !Use Runge-Kutta to compute three starting values.
      DO J=0,2
        TJ=T(J)
        YJ=Y(J)
        K1=H*F(TJ,YJ)
        K2=H*F(TJ+H/2,YJ+0.5*K1)
        K3=H*F(TJ+H/2,YJ+0.5*K2)
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
      ENDDO
      !Start of Procedure Milne-Simpson
      F1=F(T(1),Y(1))
      F2=F(T(2),Y(2))
      F3=F(T(3),Y(3))
      Pold=0
      Yold=0
      DO K=3,M-1
        Pnew=Y(K-3)+4*H*(2*F1-F2+2*F3)/3
        Pmod=Pnew+28*(Yold-Pold)/29
        T(K+1)=A+H*(K+1)
        F4=F(T(K+1),Pmod)
        Y(K+1)=Y(K-1)+H*(F2+4*F3+F4)/3
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) RETURN
        Pold=Pnew
        Yold=Y(K+1)
        F1=F2
        F2=F3
        F3=F(T(K+1),Y(K+1))
      ENDDO
      RETURN
      END

      SUBROUTINE RKF45(F,A,B,Y0,Tol,T,Y,M,MaxM,Mend)
      PARAMETER(Big=1E13)
      INTEGER I,J,K,M,Mend
      REAL A,B,T,Tol,Y,Y0
      REAL Br,H,Err,Hmin,Hmax,K1,K2,K3,K4,K5,K6,TJ,YJ
      REAL A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5,S
      REAL A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5
      REAL M1,M3,M4,M5,M6,Y1,Y2,Y3,Y4,Y5,Y6,Ygood,Ynew
      DIMENSION T(1:1000),Y(1:1000)
      EXTERNAL F
      DATA A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5,
     A     A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5
     B/ 0.25,   0.25,   0.375,   0.09375,  0.28125,   0.9230769231,
     C  0.8793809741,  -3.277196177,           3.320892126,    1.0,
     D  2.032407407,   -8.0,   7.173489279,   -0.2058966862,   0.5,
     E -0.2962962963,   2.0,  -1.381676413,    0.4529727096,  -0.275,
     F  0.002777777778,    -0.02994152047,    -0.02919989367,   0.02,
     G  0.03636363636, 0.1157407407, 0.5489278752, 0.535331384, -0.2/
      H=(B-A)/M
      Hmin=H/1024
      Hmax=H*64
      T(0)=A
      Y(0)=Y0
      T(0)=A
      J=0
      TJ=A
      Br=B-0.00001*ABS(B)
      WHILE (T(J).LT.B)
        IF ((T(J)+H).GT.Br) H=B-T(J)
        TJ=T(J)
        YJ=Y(J)
        Y1=YJ
        K1=H*F(TJ,Y1)
        Y2=YJ+B2*K1;                         IF Big.LT.ABS(Y2) RETURN
        K2=H*F(TJ+A2*H,Y2)
        Y3=YJ+B3*K1+C3*K2;                   IF Big.LT.ABS(Y3) RETURN
        K3=H*F(TJ+A3*H,Y3)
        Y4=YJ+B4*K1+C4*K2+D4*K3;             IF Big.LT.ABS(Y4) RETURN
        K4=H*F(TJ+A4*H,Y4)
        Y5=YJ+B5*K1+C5*K2+D5*K3+E5*K4;       IF Big.LT.ABS(Y5) RETURN
        K5=H*F(TJ+A5*H,Y5)
        Y6=YJ+B6*K1+C6*K2+D6*K3+E6*K4+F6*K5; IF Big.LT.ABS(Y6) RETURN
        K6=H*F(TJ+A6*H,Y6)
        Err=ABS(R1*K1+R3*K3+R4*K4+R5*K5+R6*K6)
        Ynew=YJ+N1*K1+N3*K3+N4*K4+N5*K5
        Err=ABS(Err)
        IF (Err.LT.Tol.OR.H.LE.2*Hmin) THEN
          Y(J+1)=Ynew
          IF ((TJ+H).GT.Br) THEN
              T(J+1)=B
          ELSE
              T(J+1)=TJ+H
          ENDIF
          J=J+1
          TJ=T(J)
        ENDIF
        IF (Err.EQ.0) THEN
          S=0
        ELSE
          S=0.84*(Tol*H/Err)**0.25
        ENDIF
        IF (S.LT.0.75.AND.H.GT.(2*Hmin)) H=H/2
        IF (S.GT.1.50.AND.(2*H).LT.Hmax) H=H*2
        IF (Big.LT.ABS(Y(J)).OR.MaxM.EQ.J) RETURN
        Mend=J
        IF (B.GT.T(J)) THEN
          M=J+1
        ELSE
          M=J
        ENDIF
        WRITE(9,1000) T(J),Y(J)
      REPEAT
      RETURN
1000  FORMAT(5X,2F18.9)
      END

      SUBROUTINE XHAMMING(F,A,B,Y0,M,T,Y,Mend)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      REAL A,B,Cnew,Cold,H,F1,F2,F3,F4,Pmod,Pnew,Pold,T,TK,Y,YK,Y0
      REAL K1,K2,K3,K4,TJ,YJ
      DIMENSION T(1:1000),Y(1:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
C     Use Runge-Kutta to compute three starting values.
      DO 10 J=0,2
        TJ=T(J)
        YJ=Y(J)
        K1=H*F(TJ,YJ)
        K2=H*F(TJ+H/2,YJ+0.5*K1)
        K3=H*F(TJ+H/2,YJ+0.5*K2)
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
10    CONTINUE
C     Start of Procedure Hamming
      F1=F(T(1),Y(1))
      F2=F(T(2),Y(2))
      F3=F(T(3),Y(3))
      Pold=0
      Cold=0
      DO 20 K=3,M-1
        Pnew=Y(K-3)+4*H*(2*F1-F2+2*F3)/3
        Pmod=Pnew+112*(Cold-Pold)/121
        T(K+1)=A+H*(K+1)
        F4=F(T(K+1),Pmod)
        Cnew=(9*Y(K)-Y(K-2)+3*H*(-F2+2*F3+F4))/8
        Y(K+1)=Cnew+9*(Pnew-Cnew)/121
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) RETURN
        Pold=Pnew
        Cold=Cnew
        F1=F2
        F2=F3
        F3=F(T(K+1),Y(K+1))
20    CONTINUE
      RETURN
      END

      SUBROUTINE XABM(F,A,B,Y0,M,T,Y,Mend)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      INTEGER J,K,M
      REAL A,B,H,H2,F0,F1,F2,F3,F4,P,T,TK,Y,YK,Y0
      REAL K1,K2,K3,K4,TJ,YJ
      DIMENSION T(1:1000),Y(1:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
C     Use Runge-Kutta to compute three starting values.
      DO 10 J=0,2
        TJ=T(J)
        YJ=Y(J)
        K1=H*F(TJ,YJ)
        K2=H*F(TJ+H/2,YJ+0.5*K1)
        K3=H*F(TJ+H/2,YJ+0.5*K2)
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
10    CONTINUE
C     Start of Procedure Adams-Bashforth-Moulton
      F0=F(T(0),Y(0))
      F1=F(T(1),Y(1))
      F2=F(T(2),Y(2))
      F3=F(T(3),Y(3))
      H2=H/24
      DO 20 K=3,M-1
        TK=T(K)
        YK=Y(K)
        P=YK+H2*(-9*F0+37*F1-59*F2+55*F3)
        T(K+1)=A+H*(K+1)
        F4=F(T(K+1),P)
        Y(K+1)=YK+H2*(F1-5*F2+19*F3+9*F4)
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) RETURN
        F0=F1
        F1=F2
        F2=F3
        F3=F(T(K+1),Y(K+1))
20    CONTINUE
      RETURN
      END

      SUBROUTINE XMILNE(F,A,B,Y0,M,T,Y,Mend)
C     This subroutine uses labeled DO loop(s).
      PARAMETER(Big=1E13)
      INTEGER J,K,M
      REAL A,B,H,F1,F2,F3,F4,Pmod,Pnew,Pold,T,TK,Y,YK,Yold,Y0
      REAL K1,K2,K3,K4,TJ,YJ
      DIMENSION T(1:1000),Y(1:1000)
      EXTERNAL F
      H=(B-A)/M
      T(0)=A
      Y(0)=Y0
C     Use Runge-Kutta to compute three starting values.
      DO 10 J=0,2
        TJ=T(J)
        YJ=Y(J)
        K1=H*F(TJ,YJ)
        K2=H*F(TJ+H/2,YJ+0.5*K1)
        K3=H*F(TJ+H/2,YJ+0.5*K2)
        K4=H*F(TJ+H,YJ+K3)
        Y(J+1)=YJ+(K1+2*K2+2*K3+K4)/6
        T(J+1)=A+H*(J+1)
10    CONTINUE
C     Start of Procedure Milne-Simpson
      F1=F(T(1),Y(1))
      F2=F(T(2),Y(2))
      F3=F(T(3),Y(3))
      Pold=0
      Yold=0
      DO 20 K=3,M-1
        Pnew=Y(K-3)+4*H*(2*F1-F2+2*F3)/3
        Pmod=Pnew+28*(Yold-Pold)/29
        T(K+1)=A+H*(K+1)
        F4=F(T(K+1),Pmod)
        Y(K+1)=Y(K-1)+H*(F2+4*F3+F4)/3
        Mend=K+1
        IF (Big.LT.ABS(Y(K+1))) RETURN
        Pold=Pnew
        Yold=Y(K+1)
        F1=F2
        F2=F3
        F3=F(T(K+1),Y(K+1))
20    CONTINUE
      RETURN
      END

      SUBROUTINE XRKF45(F,A,B,Y0,Tol,T,Y,M,MaxM,Mend)
      PARAMETER(Big=1E13)
      INTEGER I,J,K,M,Mend
      REAL A,B,T,Tol,Y,Y0
      REAL Br,H,Err,Hmin,Hmax,K1,K2,K3,K4,K5,K6,TJ,YJ
      REAL A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5,S
      REAL A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5
      REAL M1,M3,M4,M5,M6,Y1,Y2,Y3,Y4,Y5,Y6,Ygood,Ynew
      DIMENSION T(1:1000),Y(1:1000)
      EXTERNAL F
      DATA A2,B2,A3,B3,C3,A4,B4,C4,D4,A5,B5,C5,D5,E5,
     A     A6,B6,C6,D6,E6,F6,R1,R3,R4,R5,R6,N1,N3,N4,N5
     B/ 0.25,   0.25,   0.375,   0.09375,  0.28125,   0.9230769231,
     C  0.8793809741,  -3.277196177,           3.320892126,    1.0,
     D  2.032407407,   -8.0,   7.173489279,   -0.2058966862,   0.5,
     E -0.2962962963,   2.0,  -1.381676413,    0.4529727096,  -0.275,
     F  0.002777777778,    -0.02994152047,    -0.02919989367,   0.02,
     G  0.03636363636, 0.1157407407, 0.5489278752, 0.535331384, -0.2/
      H=(B-A)/M
      Hmin=H/1024
      Hmax=H*64
      T(0)=A
      Y(0)=Y0
      T(0)=A
      J=0
      TJ=A
      Br=B-0.00001*ABS(B)
10    IF (T(J).LT.B) THEN
        IF ((T(J)+H).GT.Br) H=B-T(J)
        TJ=T(J)
        YJ=Y(J)
        Y1=YJ
        K1=H*F(TJ,Y1)
        Y2=YJ+B2*K1
          IF (Big.LT.ABS(Y2)) RETURN
        K2=H*F(TJ+A2*H,Y2)
        Y3=YJ+B3*K1+C3*K2
          IF (Big.LT.ABS(Y3)) RETURN
        K3=H*F(TJ+A3*H,Y3)
        Y4=YJ+B4*K1+C4*K2+D4*K3
          IF (Big.LT.ABS(Y4)) RETURN
        K4=H*F(TJ+A4*H,Y4)
        Y5=YJ+B5*K1+C5*K2+D5*K3+E5*K4
          IF (Big.LT.ABS(Y5)) RETURN
        K5=H*F(TJ+A5*H,Y5)
        Y6=YJ+B6*K1+C6*K2+D6*K3+E6*K4+F6*K5
          IF (Big.LT.ABS(Y6)) RETURN
        K6=H*F(TJ+A6*H,Y6)
        Err=ABS(R1*K1+R3*K3+R4*K4+R5*K5+R6*K6)
        Ynew=YJ+N1*K1+N3*K3+N4*K4+N5*K5
        Err=ABS(Err)
        IF (Err.LT.Tol.OR.H.LE.2*Hmin) THEN
          Y(J+1)=Ynew
          IF ((TJ+H).GT.Br) THEN
              T(J+1)=B
          ELSE
              T(J+1)=TJ+H
          ENDIF
          J=J+1
          TJ=T(J)
        ENDIF
        IF (Err.EQ.0) THEN
          S=0
        ELSE
          S=0.84*(Tol*H/Err)**0.25
        ENDIF
        IF (S.LT.0.75.AND.H.GT.(2*Hmin)) THEN
          H=H/2
        ENDIF
        IF (S.GT.1.50.AND.(2*H).LT.Hmax) THEN
          H=H*2
        ENDIF
        IF (Big.LT.ABS(Y(J)).OR.MaxM.EQ.J) RETURN
        Mend=J
        IF (B.GT.T(J)) THEN
          M=J+1
        ELSE
          M=J
        ENDIF
        WRITE(9,1000) T(J),Y(J)
        GOTO 10
      ENDIF
      RETURN
1000  FORMAT(5X,2F18.9)
      END

      SUBROUTINE MESSAGE(Meth,Order,State)
      INTEGER K,Meth,Order,State
      REAL Resp
      WRITE(9,*)'          Solution of the differential equation  Y` = F
     +(T,Y)'
      WRITE(9,*)' '
      WRITE(9,*)'          with the initial condition  Y(A) = Y .'
      WRITE(9,*)'                                              0'
      WRITE(9,*)' '
      WRITE(9,*)'          A numerical approximation is computed over [A
     +,B].'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'          Choose the method of approximation:'
      WRITE(9,*)' '
      WRITE(9,*)'    < 1 > Adams-Bashforth-Moulton method'
      WRITE(9,*)' '
      WRITE(9,*)'    < 2 > Milne-Simpson`s method'
      WRITE(9,*)' '
      WRITE(9,*)'    < 3 > Hamming`s method'
      WRITE(9,*)' '
      WRITE(9,*)'    < 4 > Taylor`s method of order N=3'
      WRITE(9,*)' '
      WRITE(9,*)'    < 5 > Taylor`s method of order N=4'
      WRITE(9,*)' '
      WRITE(9,*)'    < 6 > Runge-Kutta-Fehlberg'
      WRITE(9,*)' '
      WRITE(9,*)'          SELECT < 1 - 6 > ? '
      Resp=Meth
      READ(9,*) Resp
      Meth=INT(Resp)
      IF ((Meth.LT.1.OR.Meth.GT.6).AND.State.NE.0) THEN
        Meth=1
      ENDIF
      IF (Meth.EQ.4) THEN
        Order=3
      ENDIF
      IF (Meth.EQ.5) THEN
        Order=4
      ENDIF
      RETURN
      END

      SUBROUTINE INPUT(Meth)
      INTEGER Meth
      REAL Valu
      CHARACTER Ans*1
      DO 10 I=1,15
        WRITE(9,*)' '
10    CONTINUE
      WRITE(9,*)'     You chose '
      WRITE(9,*)' '
      IF (Meth.EQ.1) THEN
        WRITE(9,*)'     Adams-Bashforth-Moulton method to solve
     + the D.E.'
      ELSEIF (Meth.EQ.2) THEN
        WRITE(9,*)'     the Milne-Simpson method to solve
     + the D.E.'
      ELSEIF (Meth.EQ.3) THEN
        WRITE(9,*)'     Hamming`s method to solve the D.E.'
      ELSEIF (Meth.EQ.4) THEN
        WRITE(9,*)'     Taylor`s method of order N=3 to solve
     + the D.E.'
      ELSEIF (Meth.EQ.5) THEN
        WRITE(9,*)'     Taylor`s method of order N=4 to solve
     + the D.E.'
      ELSEIF (Meth.EQ.6) THEN
        WRITE(9,*)'     the Runge-Kutta-Fehlberg method to solve
     + the D.E.'
      ENDIF
      WRITE(9,*)' '
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)'     with the initial condition Y(A) = Y0.'
      WRITE(9,*)' '
      WRITE(9,*)'     A numerical approximation is computed over [A,B].'
      WRITE(9,*)' '
      WRITE(9,*)'     You must enter the endpoints for the interval,'
      WRITE(9,*)' '
      WRITE(9,*)'     the initial condition Y0, and the number of steps
     + M.'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'         Press the <ENTER> key. '
      READ(9,'(A)') Ans
      RETURN
      END

      SUBROUTINE EPOINTS(A,B,Y0,M,Tol,State)
      INTEGER I,M,Stat,State
      REAL A,B,Tol,Y0,Valu
C     STATUS = (Change,Enter,Done)
      CHARACTER Ans*1
      Stat=1
      IF (State.EQ.0) THEN
        Stat=0
      ENDIF
10    CONTINUE
      DO 20 I=1,15
        WRITE(9,*)' '
20    CONTINUE
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)' '
      IF (Stat.EQ.1) THEN
        WRITE(9,*)'     ENTER  the  left  endpoint  A = '
        READ(9,*) A
        WRITE(9,*)' '
        WRITE(9,*)'     ENTER  the  right endpoint  B = '
        READ(9,*) B
        WRITE(9,*)' '
        WRITE(9,*)'     ENTER initial condition  Y(A) = '
        READ(9,*) Y0
        WRITE(9,*)' '
        WRITE(9,*)'     ENTER the number of steps   M = '
        Valu=M
        READ(9,*) Valu
        M=INT(Valu)
        IF (M.LT.1) THEN
          M=1
        ENDIF
        IF (M.GT.1000) THEN
          M=1000
        ENDIF
        WRITE(9,*)'     ENTER   the   tolerance   Tol = '
        Valu=Tol
        READ(9,*) Valu
        Tol=Valu
      ELSE
        WRITE(9,*)'     The left  endpoint is    A =',A
        WRITE(9,*)' '
        WRITE(9,*)'     The right endpoint is    B =',B
        WRITE(9,*)' '
        WRITE(9,*)'     Initial condition is  Y(A) =',Y0
        WRITE(9,*)' '
        WRITE(9,*)'     The number of steps is   M = ',M
        WRITE(9,*)' '
        WRITE(9,*)'     The   tolerance  is    Tol = ',Tol
      ENDIF
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'     Do you want to make a change ?  <Y/N>  '
      READ(9,'(A)') Ans
      IF (Ans.EQ.'Y'.OR.Ans.EQ.'y') THEN
          Stat=0
          DO 30 I=1,15
            WRITE(9,*)' '
30        CONTINUE
          WRITE(9,*)'          '
          CALL PRINTFUN
          WRITE(9,*)' '
          WRITE(9,*)' '
          WRITE(9,*)'     The current left  endpoint is A =',A
            WRITE(9,*)'     ENTER  the NEW left  endpoint A = '
          READ(9,*) A
          WRITE(9,*)' '
          WRITE(9,*)'     The current right endpoint is B =',B
            WRITE(9,*)'     ENTER  the NEW right endpoint B = '
          READ(9,*) B
          WRITE(9,*)' '
          WRITE(9,*)'        The current   I. C. is  Y(A) =',Y0
            WRITE(9,*)'        ENTER the NEW    I. C.  Y(A) = '
          READ(9,*) Y0
          WRITE(9,*)' '
          WRITE(9,*)'        The current value of M is  M = ',M
            WRITE(9,*)'        ENTER  the  NEW  value of  M = '
          Valu=M
          READ(9,*) Valu
          M=INT(Valu)
          IF (M.LT.1) THEN
            M=1
          ENDIF
          IF (M.GT.1000) THEN
            M=1000
          ENDIF
          WRITE(9,*)'    The current value of Tol is  Tol = ',Tol
            WRITE(9,*)'        ENTER  the NEW value of  Tol = '
          READ(9,*) Valu
          Tol=Valu
      ELSE
        Stat=2
      ENDIF
      IF (Stat.EQ.0.OR.Stat.EQ.1) GOTO 10
      RETURN
      END

      SUBROUTINE RESULTS(T,Y,M,Mend,Meth)
      INTEGER I,K,M,Mend,Meth
      REAL T,Y
      DIMENSION T(1:1000),Y(1:1000)
      DO 10 I=1,10
        WRITE(9,*)' '
10    CONTINUE
      IF (Meth.EQ.1) THEN
        WRITE(9,*)'Adams-Bashforth-Moulton method was used to solve the
     & D.E.'
      ELSEIF (Meth.EQ.2) THEN
        WRITE(9,*)'The Milne-Simpson method was used to solve the D.E.'
      ELSEIF (Meth.EQ.3) THEN
        WRITE(9,*)'Hamming`s method was used to solve the D.E.'
      ELSEIF (Meth.EQ.4) THEN
        WRITE(9,*)'Taylor`s method of order N=3 was used to solve the
     & D.E.'
      ELSEIF (Meth.EQ.5) THEN
        WRITE(9,*)'Taylor`s method of order N=4 was used to solve the
     & D.E.'
      ELSEIF (Meth.EQ.6) THEN
        WRITE(9,*)'The Runge-Kutta-Fehlberg method to solve the D.E.'
      ENDIF
      WRITE(9,*)' '
      CALL PRINTFUN
      WRITE(9,*)' '
      WRITE(9,*)'with  Y(',T(0),'  ) =',Y(0)
      WRITE(9,*)' '
      WRITE(9,*)'   K            T(K)           Y(K)'
      WRITE(9,*)' ---------------------------------------'
      DO 20 K=0,Mend
        WRITE(9,999) K, T(K), Y(K)
20    CONTINUE
      IF (Mend.LT.M) THEN
        WRITE(9,*)' '
        WRITE(9,*)'The solution points are approaching a pole.'
      ENDIF
999   FORMAT(I4,4X,4F15.7)
      PAUSE
      RETURN
      END

      PROGRAM BDDVALUE
C     NUMERICAL METHODS: FORTRAN Programs, (c) John H. Mathews 1994
C     To accompany the text:
C     NUMERICAL METHODS for Mathematics, Science and Engineering, 2nd Ed, 1992
C     Prentice Hall, Englewood Cliffs, New Jersey, 07632, U.S.A.
C     This free software is complements of the author.
C
C     Algorithm 9.9 (Linear Shooting Method).
C     Section 9.8, Boundary Value Problems, Page 488
C
C     Algorithm 9.10 (Finite-Difference Method).
C     Section 9.9, Finite-Difference Method, Page 496
C
      PARAMETER(GNmax = 630,MaxM = 101)
      CHARACTER Ans*1,Ach*1,Mess*80
      REAL D,T,X,X1,X2
      DIMENSION D(1:4),T(0:MaxM),X(0:MaxM),X1(0:MaxM),X2(0:MaxM)
      INTEGER FunM,FunT, GNpts, Inum, M, Mend, Meth, Order, Sub
      REAL A, B, Alpha, Beta
      INTEGER State,DoMo,Changes,Done,Working,Go,Stop
      EXTERNAL F1,F2,P,Q,R
      FunM = 9
      Meth = 1
      CALL MESSAGE(Meth)
      CALL INPUT(Meth,FunT,FunM)
      CALL EPOINTS(A,B,Alpha,Beta,M)
      IF Meth.EQ.1 THEN
        CALL LinSho(A,B,Alpha,Beta,M,T,X,FunT)
      ENDIF
      IF Meth.EQ.2 THEN
        CALL FinDif(A,B,Alpha,Beta,M,T,X,FunT)
      ENDIF
      CALL RESULTS(FunT,Meth,T,X,M,Mend)
      WRITE(9,*)' '
      READ(9,*) Ans
      STOP
      END !End of the Main Program

      REAL FUNCTION P(T,FunT)
      INTEGER FunT
      REAL T
      IF (FunT.EQ.1) THEN
        P = 2.0*T/(1.0 + T*T)
      ELSEIF (FunT.EQ.2) THEN
        P = -2.0/T
      ELSEIF (FunT.EQ.3) THEN
        P = -5.0
      ELSEIF (FunT.EQ.4) THEN
        P = -2.0
      ELSEIF (FunT.EQ.5) THEN
        P = -1.0/T
      ELSEIF (FunT.EQ.6) THEN
        P = -2.0/T
      ELSEIF (FunT.EQ.7) THEN
        P = -1.0/T
      ELSEIF (FunT.EQ.8) THEN
        P = 1.0/T
      ELSEIF (FunT.EQ.9) THEN
        P = -1.0/T
      ENDIF
      RETURN
      END

      REAL FUNCTION Q(T,FunT)
      INTEGER FunT
      REAL T
      IF (FunT.EQ.1) THEN
        Q = -2.0/(1.0 + T*T)
      ELSEIF (FunT.EQ.2) THEN
        Q = 2.0/(T*T)
      ELSEIF (FunT.EQ.3) THEN
        Q = -6.0
      ELSEIF (FunT.EQ.4) THEN
        Q = -2.0
      ELSEIF (FunT.EQ.5) THEN
        Q = -(1.0 - 1.0/(4.0*T*T))
      ELSEIF (FunT.EQ.6) THEN
        Q = 2.0/(T*T)
      ELSEIF (FunT.EQ.7) THEN
        Q = -(1.0 - 1.0/(4.0*T*T))
      ELSEIF (FunT.EQ.8) THEN
        Q = -1.0/(T*T)
      ELSEIF (FunT.EQ.9) THEN
        Q = -16.0/(T*T)
      ENDIF
      RETURN
      END

      REAL FUNCTION R(T,FunT)
      INTEGER FunT
      REAL T
      IF (FunT.EQ.1) THEN
        R = 1.0
      ELSEIF (FunT.EQ.2) THEN
        R = 10.0*COS(LOG(T))/(T*T)
      ELSEIF (FunT.EQ.3) THEN
        R = T*EXP(-2.0*T) + 3.9*COS(3.0*T)
      ELSEIF (FunT.EQ.4) THEN
        R = EXP(-T) + SIN(2.0*T)
      ELSEIF (FunT.EQ.5) THEN
        R = 0.0
      ELSEIF (FunT.EQ.6) THEN
        R = SIN(T)/(T*T)
      ELSEIF (FunT.EQ.7) THEN
        R = SQRT(T)*COS(T)
      ELSEIF (FunT.EQ.8) THEN
        R = 1.0
      ELSEIF (FunT.EQ.9) THEN
        R = 1.0/(T*T)
      ENDIF
      RETURN
      END

      REAL FUNCTION F1(T,X,Y,FunT)
      INTEGER FunT
      REAL T,X,Y
      IF (FunT.EQ.1) THEN
        F1 = 2.0*T*Y/(1.0 + T*T) - 2.0*X/(1.0 + T*T) + 1.0
      ELSEIF (FunT.EQ.2) THEN
        F1 = -2.0/T*Y + 2.0/(T*T)*X + 10.0*COS(LOG(T))/(T*T)
      ELSEIF (FunT.EQ.3) THEN
        F1 = -5.0*Y - 6.0*X + T*EXP(-2.0*T) + 3.9*COS(3.0*T)
      ELSEIF (FunT.EQ.4) THEN
        F1 = -2.0*Y - 2.0*X + EXP(-T) + SIN(2.0*T)
      ELSEIF (FunT.EQ.5) THEN
        F1 = -1.0/T*Y - (1.0 - 1.0/(4.0*T*T))*X
      ELSEIF (FunT.EQ.6) THEN
        F1 = -2.0/T*Y + 2.0/(T*T)*X + SIN(T)/(T*T)
      ELSEIF (FunT.EQ.7) THEN
        F1 = -1.0/T*Y-(1.0 - 1.0/(4.0*T*T))*X+SQRT(T)*COS(T)
      ELSEIF (FunT.EQ.8) THEN
        F1 = 1.0/T*Y - 1.0/(T*T)*X + 1.0
      ELSEIF (FunT.EQ.9) THEN
        F1 = -1.0/T*Y - 16.0/(T*T)*X + 1.0/(T*T)
      ENDIF
      RETURN
      END

      REAL FUNCTION F2(T,X,Y,FunT)
      INTEGER FunT
      REAL T,X,Y
      IF (FunT.EQ.1) THEN
        F2 = 2.0*T*Y/(1.0 + T*T) - 2.0*X/(1.0 + T*T)
      ELSEIF (FunT.EQ.2) THEN
        F2 = -2.0/T*Y + 2.0/(T*T)*X
      ELSEIF (FunT.EQ.3) THEN
        F2 = -5.0*Y - 6.0*X
      ELSEIF (FunT.EQ.4) THEN
        F2 = -2.0*Y - 2.0*X
      ELSEIF (FunT.EQ.5) THEN
        F2 = -1.0/T*Y - (1.0 - 1.0/(4.0*T*T))*X
      ELSEIF (FunT.EQ.6) THEN
        F2 = -2.0/T*Y + 2.0/(T*T)*X
      ELSEIF (FunT.EQ.7) THEN
        F2 = -1.0/T*Y - (1.0 - 1.0/(4.0*T*T))*X
      ELSEIF (FunT.EQ.8) THEN
        F2 = 1.0/T*Y - 1.0/(T*T)*X
      ELSEIF (FunT.EQ.9) THEN
        F2 = -1.0/T*Y - 16/(T*T)*X
      ENDIF
      RETURN
      END

      SUBROUTINE PRTFUN(FunT)
      INTEGER FunT
      IF (FunT.EQ.1) THEN
        WRITE(9,*)'X`` =  2*T*X`/(1 + T*T) - 2*X/(1 + T*T) + 1'
      ELSEIF (FunT.EQ.2) THEN
        WRITE(9,*)'X`` =  -2/T*X` + 2/(T*T)*X + 10*COS(LOG(T))/(T*T)'
      ELSEIF (FunT.EQ.3) THEN
        WRITE(9,*)'X`` =  -5*X` - 6*X + T*EXP(-2*T) + 3.9*COS(3*T)'
      ELSEIF (FunT.EQ.4) THEN
        WRITE(9,*)'X`` =  -2*X` - 2*X + EXP(-T) + SIN(2*T)'
      ELSEIF (FunT.EQ.5) THEN
        WRITE(9,*)'X`` =  -1/T*X` - (1 - 1/(4*T*T))*X'
      ELSEIF (FunT.EQ.6) THEN
        WRITE(9,*)'X`` =  -2/T*X` + 2/(T*T)*X + SIN(T)/(T*T)'
      ELSEIF (FunT.EQ.7) THEN
        WRITE(9,*)'X`` =  -1/T*X` - (1 - 1/(4*T*T))*X + SQRT(T)*COS(T)'
      ELSEIF (FunT.EQ.8) THEN
        WRITE(9,*)'X`` =  1/T*X` - 1/(T*T)*X + 1'
      ELSEIF (FunT.EQ.9) THEN
        WRITE(9,*)'X`` =   -1/T*X` - 16/(T*T)*X + 1/(T*T)'
      ENDIF
      RETURN
      END

      !Start of the subroutine for the Linear Shooting method.
      SUBROUTINE RKbdF(F,A,B,Alpha,Beta,M,T,X,FunT)
      PARAMETER(MaxM = 101)
      INTEGER J,M,FunT
      REAL A,B,Alpha,Beta,T,X,Y
      REAL H,S1,S2 S3,S4,R1,R2,R3,R4,Tj,Xj,Yj
      DIMENSION T(0:MaxM),X(0:MaxM),Y(0:MaxM)
      EXTERNAL F
      H = (B - A)/M
      T(0) = A
      X(0) = Alpha
      Y(0) = Beta
      DO J=0,(M-1)
        Tj = T(J)
        Xj = X(J)
        Yj = Y(J)
        S1 = H*Yj
        R1 = H*F(Tj, Xj, Yj,FunT)
        S2 = H*(Yj + R1/2.0)
        R2 = H*F(Tj + H/2.0, Xj+S1/2.0, Yj+R1/2.0,FunT)
        S3 = H*(Yj + R2/2.0)
        R3 = H*F(Tj + H/2.0, Xj+S2/2.0, Yj+R2/2.0,FunT)
        S4 = H*(Yj + R3)
        R4 = H*F(Tj + H, Xj + S3, Yj + R3,FunT)
        X(J+1) = Xj + (S1 + 2.0*S2 + 2.0*S3 + S4)/6.0
        Y(J+1) = Yj + (R1 + 2.0*R2 + 2.0*R3 + R4)/6.0
        T(J+1) = A + H*(J+1)
      ENDDO
      END

      !Start of the main program for the Linear Shooting method.
      SUBROUTINE LinSho(A,B,Alpha,Beta,M,T,X,FunT)
      PARAMETER(MaxM = 101)
      INTEGER J,M,FunT
      REAL A,B,Alpha,Beta
      REAL T,X,X1,X2
      DIMENSION T(0:MaxM),X(0:MaxM),X1(0:MaxM),X2(0:MaxM)
      EXTERNAL F1,F2
      CALL RKbdF(F1,A,B,Alpha,0.0,M,T,X1,FunT)
      CALL RKbdF(F2,A,B, 0.0 ,1.0,M,T,X2,FunT)
      DO J=0,M
        X(J) = X1(J) + (Beta - X1(M))*X2(J)/X2(M)
      ENDDO
      END   !The main program for the Linear Shooting method.

      !Start of the subroutines for the Finite Difference method.
      SUBROUTINE CoeMat(A,B,Alpha,Beta,N,Va,Vb,Vc,Vd,Vt,FunT)
      PARAMETER(MaxM = 101)
      INTEGER J,N,FunT
      REAL A,B,Alpha,Beta
      REAL Va, Vb, Vc, Vd, Vt
      DIMENSION Va(0:MaxM),Vb(0:MaxM),Vc(0:MaxM),Vd(0:MaxM),Vt(0:MaxM)
      EXTERNAL P,Q,R
      H = (B - A)/N
      DO J=1,(N-1)
        Vt(J) = A + H*J
      ENDDO
      DO J=1,(N-1)
        Vb(J) = -H*H*R(Vt(J),FunT)
      ENDDO
      Vb(1) = Vb(1) + (1 + H/2*P(Vt(1),FunT))*Alpha
      Vb(N-1) = Vb(N-1) + (1 - H/2*P(Vt(N-1),FunT))*Beta
      DO J=1,(N-1)
        Vd(J) = 2 + H*H*Q(Vt(J),FunT)
      ENDDO
      DO J=1,(N-2)
        Va(J) = -1 - H/2*P(Vt(J+1),FunT)
      ENDDO
      DO J=1,(N-2)
        Vc(J) = -1 + H/2*P(Vt(J),FunT)
      ENDDO
      END

      SUBROUTINE TriMat(Va,Vb,Vc,Vd,N,X)
      PARAMETER(MaxM = 101)
      INTEGER K,N
      REAL Va,Vb,Vc,Vd,X
      DIMENSION Va(0:MaxM),Vb(0:MaxM),Vc(0:MaxM),Vd(0:MaxM),X(0:MaxM)
      REAL TEM
      DO K=2,(N-1)
        TEM = Va(K-1)/Vd(K-1)
        Vd(K) = Vd(K) - TEM*Vc(K-1)
        Vb(K) = Vb(K) - TEM*Vb(K-1)
      ENDDO
      X(N-1) = Vb(N-1)/Vd(N-1)
      DO K=(N-2),1,-1
        X(K) = (Vb(K) - Vc(K)*X(K+1))/Vd(K)
      ENDDO
      END

      !Start of the main program for the Finite Difference method.
      SUBROUTINE FinDif(A,B,Alpha,Beta,N,T,X,FunT)
      PARAMETER(MaxM = 101)
      INTEGER N,FunT
      REAL A,B,Alpha,Beta,H
      REAL T, X
      DIMENSION T(0:MaxM), X(0:MaxM)
      REAL Va, Vb, Vc, Vd
      DIMENSION Va(0:MaxM), Vb(0:MaxM), Vc(0:MaxM), Vd(0:MaxM)
      !Remark: Some of the zero components are NOT used.
      EXTERNAL P,Q,R
      CALL CoeMat(A,B,Alpha,Beta,N,Va,Vb,Vc,Vd,T,FunT)
      CALL TriMat(Va, Vb, Vc, Vd, N, X)
      T(0) = A
      T(N) = B
      X(0) = Alpha
      X(N) = Beta
      END   !The main program for the Finite Difference method.

      SUBROUTINE MESSAGE(Meth)
      INTEGER Meth
      CHARACTER Mess*80
      WRITE(9,*)'       BOUNDARY VALUE DIFFERENTIAL EQUATIONS'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'Solution of the differential equation'
      WRITE(9,*)' '
      WRITE(9,*)'x`` = p(t)x ` ( t ) + q ( t ) x ( t ) + r ( t )'
      WRITE(9,*)' '
      WRITE(9,*)'with x(a)=Alpha and x(b)=Beta computed over [a,b].'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'       Choose the method of approximation:'
      WRITE(9,*)' '
      WRITE(9,*)'       < 1 >  Linear shooting method'
      WRITE(9,*)' '
      WRITE(9,*)'       < 2 >  Finite difference method'
      WRITE(9,*)' '
      Mess = '       SELECT < 1 - 2 > ?  '
      Meth = 1
      WRITE(9,*) Mess
      READ(9,*)  Meth
      IF (Meth .LE. 1) .AND. (State .NE. Changes) THEN
        Meth = 1
      ENDIF
      IF (Meth .GT. 2) .AND. (State .NE. Changes) THEN
        Meth = 2
      ENDIF
      END

      SUBROUTINE INPUT(Meth,FunT,FunM)
      INTEGER FunT,FunM,K
      CHARACTER Ans*1,Mess*80
      WRITE(9,*)' '
      IF Meth.EQ.1 THEN
        WRITE(9,*)'You chose the linear shooting method'
      ENDIF
      IF Meth.EQ.2 THEN
        WRITE(9,*)'You chose the finite difference method'
      ENDIF
      WRITE(9,*)'to solve Y` = F(T,Y).'
      WRITE(9,*)' '
      WRITE(9,*)'Choose your D.E.:'
      WRITE(9,*)' '
      DO K=1,FunM
        WRITE(9,*)'<', K,'>'
        CALL PRTFUN(K)
        WRITE(9,*)' '
      ENDDO
      Mess = '     SELECT < 1 - 9 > ?  '
      WRITE(9,*) Mess
      READ(9,*) FunT
      IF (FunT.LT.1) THEN
        FunT = 1
      ENDIF
      IF (FunT.GT.FunM) THEN
        FunT = FunM
      ENDIF
      WRITE(9,*)' '
      IF Meth.EQ.1 THEN
        WRITE(9,*)'You chose the linear shooting method'
      ENDIF
      IF Meth.EQ.2 THEN
        WRITE(9,*)'You chose the finite difference method'
      ENDIF
      WRITE(9,*)'to solve the D.E.'
      WRITE(9,*)' '
      WRITE(9,*)'                     '
      CALL PRTFUN(FunT)
      WRITE(9,*)' '
      WRITE(9,*)'with the initial conditions'
      WRITE(9,*)' '
      WRITE(9,*)'X(A) = Alpha  and  X(B) = Beta.'
      WRITE(9,*)' '
      WRITE(9,*)'A numerical approximation is computed over [A,B].'
      WRITE(9,*)' '
      WRITE(9,*)'You must give the endpoints for the interval,'
      WRITE(9,*)' '
      WRITE(9,*)'the boundary values, and the number of steps.'
      WRITE(9,*)' '
      WRITE(9,*)' '
      WRITE(9,*)'     Press the <ENTER> key. '
      READ(9,*) Ans
      WRITE(9,*)' '
      end

      SUBROUTINE EPOINTS(A,B,Alpha,Beta,M)
      INTEGER M
      REAL A,B,Alpha,Beta
      CHARACTER Mess*80
      WRITE(9,*)' '
      WRITE(9,*)'                  '
      CALL PRTFUN(FunT)
      WRITE(9,*)' '
      WRITE(9,*)' '
      Mess = 'ENTER  the  left  endpoint  A = '
      WRITE(9,*) Mess
      READ(9,*)  A
      WRITE(9,*)' '
      Mess = 'ENTER  the  right endpoint  B = '
      WRITE(9,*) Mess
      READ(9,*)  B
      WRITE(9,*)' '
      Mess = 'ENTER initial condition  X(A) = '
      WRITE(9,*) Mess
      READ(9,*) Alpha
      WRITE(9,*)' '
      Mess = 'ENTER initial condition  X(B) = '
      WRITE(9,*) Mess
      READ(9,*)  Beta
      WRITE(9,*)' '
      Mess = 'ENTER the number of steps   M = '
      M = 1
      WRITE(9,*) Mess
      READ(9,*) M
      WRITE(9,*)' '
      IF (M.LT.1) M = 1
      IF (M.GT.1000) M = 1000
      END

      SUBROUTINE RESULTS(FunT,Meth,T,X,M,Mend)
      INTEGER FunT,K,M,Mend
      REAL T,X
      DIMENSION T(0:MaxM),X(0:MaxM)
      WRITE(9,*)' '
      IF Meth.EQ.1 THEN
        WRITE(9,*)' The linear shooting method'
      ENDIF
      IF Meth.EQ.2 THEN
        WRITE(9,*)' The finite difference method'
      ENDIF
      WRITE(9,*)' was used to solve the D.E.'
      WRITE(9,*)'       '
      CALL PRTFUN(FunT)
      WRITE(9,*)' '
      WRITE(9,*)'with  X(', T(0), '  ) =', X(0)
      WRITE(9,*)'and   X(', T(M), '  ) =', X(M)
      WRITE(9,*)' '
      WRITE(9,*)'    K', '         T(K)    ', '         X(K)'
      WRITE(9,*)'  --------------------------------------------'
      WRITE(9,*)' '
      DO K=0,M
        WRITE(9,*)'  ',K, '        ', T(K), '        ', X(K)
      ENDDO
      Mend = M
      IF (Mend.LT.M) THEN
        WRITE(9,*)'The solution points are approaching a pole.'
        WRITE(9,*)' '
      ENDIF
      WRITE(9,*)'                   Press the <ENTER> key.  '
      READ(9,*) Ans
      WRITE(9,*)' '
      END