C       RUNGE-KUTTA FOR SYSTEMS BY DR. J. W. ROGERS, Auburn University
C       THIS PROGRAM USES A K-K 4TH ORDER ALGORITHM TO SOLVE THE
C       SYSTEM OF DIFFERENTIAL EQUATIONS X'=F(X). THE PARAMETER
C       STATEMENT RESERVES SPACE FOR STATE VECTORS OF DIMENSION UP TO
C       10. THE CHECKING ROUTINE USES DIM=2.
C
C
        PROGRAM MAIN
        INTEGER*2 DIMLIM
        PARAMETER(DIMLIM = 10)
        INTEGER*2 DIM,I,NUM
        REAL*8 DT,PI,T,X
        EXTERNAL GETDX
        DIMENSION X(DIMLIM)
        DATA DIM/2/

        PI = 4.D0*ATAN(1.D0)
        X(1) = 1.D0
        X(2) = 1.D0
        WRITE(*,*) 'ENTER THE NUMBER OF STEPS : '
        READ(*,*) NUM
        DT = 0.1D0
        T = 0.D0
        DO 10 I=1,10
                     CALL RKUPDATE(DIM,T,X,GETDX,DT,X)
                     T = T + DT
                     WRITE(*,1) T,X(1),X(2)
10      CONTINUE
        WRITE(*,1) T,X(1),X(2)
1       FORMAT(' T=',F8.6,2F15.10)
        END
C
C	-------------------------------------------------------------------
C
C	Here is where your equations of motion go.
C
        SUBROUTINE GETDX(T,X,DX)
        INTEGER*2 DIMLIM
        PARAMETER(DIMLIM=10)
        REAL*8 DX,T,X
        DIMENSION DX(DIMLIM), X(DIMLIM)
        DX(1) = X(2)
        DX(2) = 2.71828182846**T
        RETURN
        END
C
C	-------------------------------------------------------------------
C
C	RUNGE-KUTTA 4th ORDER
C
C       DIM: DIMENSION OF THE STATE VECTOR
C       CURRX: CURRENT STATE VECTOR
C       NEWX: STATE VECTOR AFTER UPDATE
C       RK(SBSTP,*): K VALUES FOR SUBSTEP 1,2,3,4
C       GETDX: SUBROUTINE THAT COMPUTES THE DERIVATIVE OF X
C
C	The only thing you should have to modify in this routine is the passing
C	of the parameter MU (earth's gravitation constant) to this routine and
C	then sending MU on to the subroutine GETDX.
C
        SUBROUTINE RKUPDATE(DIM,T,CURRX,GETDX,DT,NEWX)
        INTEGER*2 DIMLIM
        PARAMETER(DIMLIM=5)
        INTEGER*2 DIM
        REAL*8 DT,CURRX(DIMLIM),NEWX(DIMLIM),T
        EXTERNAL GETDX
        INTEGER*2 I,SBSTP
        REAL*8 DX(DIMLIM),OFFST(4),RK(4,DIMLIM),T1,X1(DIMLIM),XX
        SAVE OFFST
        DATA OFFST/0.D0, .5D0, .5D0, 1.D0/
C       COMPUTE K VALUES FOR RUNGE-KUTTA STEP
        DO 40 SBSTP=1,4
                     IF(SBSTP .EQ. 1) THEN
                          DO 10 I=1,DIM
                                X1(I) = CURRX(I)
10                        CONTINUE
                     ELSE
                          DO 20 I=1,DIM
                                X1(I) = CURRX(I) +
     &                                  OFFST(SBSTP)*RK(SBSTP-1,I)
20                        CONTINUE
                     ENDIF
                     T1 = T + OFFST(SBSTP) * DT
C       COMPUTE X1 AND NEXT K VALUE
                     CALL GETDX(T1,X1,DX)
                     DO 30 I=1,DIM
                          RK(SBSTP,I) = DT * DX(I)
30                   CONTINUE
40      CONTINUE
C       UPDATE STATE VECTOR
        DO 50 I=1,DIM
                     XX = RK(1,I) + 2.D0*(RK(2,I)+RK(3,I))+RK(4,I)
                     NEWX(I) = CURRX(I) + XX/6.D0
50      CONTINUE
        RETURN
        END