      SUBROUTINE RK(N,T,H,A,F)
C     ******************************************************************
C     FOURTH ORDER RUNGE-KUTTA PROCEDURE FOR SOLVING DY=F(T,Y)
C     ******************************************************************
      DIMENSION A(1)
      EXTERNAL F
      NP1=N+1
      IF (H.EQ.0.0) GO TO 50
C
      HA=.5*H
      TA=T+HA
      M=N+N
      M1=M+1
C
      DO 10 K=1,N
      NK=N+K
      MK=M+K
      A(NK)=HA*A(NK)
   10 A(MK)=A(K)+A(NK)
      CALL F(TA,A(M1))
C
      DO 20 K=1,N
      NK=N+K
      MK=M+K
      A(NK)=A(NK)+H*A(MK)
   20 A(MK)=A(K)+HA*A(MK)
      CALL F(TA,A(M1))
C
      T=T+H
      DO 30 K=1,N
      NK=N+K
      MK=M+K
      A(MK)=H*A(MK)
      A(NK)=A(NK)+A(MK)
   30 A(MK)=A(K)+A(MK)
      CALL F(T,A(M1))
C
      DO 40 K=1,N
      NK=N+K
      MK=M+K
      A(K)=A(K)+(A(NK)+HA*A(MK))/3.0
   40 A(NK)=A(K)
      CALL F(T,A(NP1))
      RETURN
C
   50 DO 51 K=1,N
      NK=N+K
   51 A(NK)=A(K)
      CALL F(T,A(NP1))
      RETURN
      END
