      PROGRAM RK_STABILITY
C...
C...  PROGRAM RK_STABILITY COMPUTES THE STABILITY INTERVAL ALONG
C...  THE NEGATIVE REAL AXIS OF AN EXPLICIT RUNGE KUTTA METHOD OF
C...  ORDER N (THE PRODUCT LAMBDA*H IS PLOTTED ALONG THE NEGATIVE
C...  REAL AXIS WHERE LAMBDA IS THE EIGENVALUE OF THE FIRST ORDER
C...  ODE DY/DT = LAMBDA*Y AND H IS THE INTEGRATION STEP USED IN 
C...  THE RUNGE KUTTA METHOD).
C...
C...  REAL SELECTED VARIABLES (LAMH IS LAMBDA*H)
      REAL LAMH 
C...
C...  OPEN AN OUTPUT FILE
      OPEN(6,FILE='OUTPUT',STATUS='NEW')
C...
C...  SET THE ORDER OF THE RK METHOD
      N=2
C...
C...  PRINT A HEADING FOR THE NUMERICAL OUTPUT
         WRITE(6,4)N
4        FORMAT(' ORDER OF THE EXPLICIT RUNGE KUTTA METHOD = ',I2,//,
     +      10X,' LAMBDA*H',5X,'Y(I+1)/Y(I)',/)
C...
C...  STEP THROUGH A SERIES OF VALUES OF LAMBDA*H ALONG THE NEGATIVE
C...  REAL AXIS FROM LAMDA*H = 0 TO LAMBDA*H = -4.  AT EACH VALUE OF
C...  LAMBDA*H, COMPUTE THE AMPLIFICATION FACTOR, I.E., Y(I+1)/Y(I),
C...  FOR THE RK METHOD
C...
      STEP=0.25
      LAMH=STEP
      DO 1 J=1,17
         LAMH=LAMH-STEP
         AMP=1.0
         DO 2 I=1,N
            AMP=AMP+TERM(I,LAMH)
2        CONTINUE
C...
C...  PRINT LAMBDA*H AND THE AMPLIFICATION FACTOR FOR THE RK METHOD
         WRITE(6,3)LAMH,AMP
3        FORMAT(2X,F15.2,F15.4)
C...
C...  GO ON TO THE NEXT LAMBDA*H
1     CONTINUE
C...
C...  ALL OF THE LAMBDA*H VALUES HAVE BEEN USED
      STOP
      END
      REAL FUNCTION TERM(I,LAMH)
C...
C...  FUNCTION TERM COMPUTES THE ITH TERM IN THE TAYLOR SERIES
C...  EXPANSION OF EXP(LAMDBA*H)
C...
      REAL LAMH
      TERM=LAMH**I/FACT(I)
      RETURN
      END
      REAL FUNCTION FACT(I)
C...
C...  FUNCTION FACT COMPUTES THE ITH FACTORIAL IN THE TAYLOR SERIES
C...  EXPANSION OF LAMBDA*H
C...
      FACT=1.0
      DO 1 J=1,I
         FACT=FACT*FLOAT(J)
1     CONTINUE
      RETURN
      END
