C #RKSTB   V1A   19-JAN-73
C LAST UPDATE:
C
      SUBROUTINE RKSTB(DERIV,Y0,YP0,Y1,N,X,H,IX,WK)
      DIMENSION Y0(N), YP0(N), Y1(N), WK(N,3)
C
C *PURPOSE.
C 	TO MAKE A SINGLE STEP IN INTEGRATING A SET OF ORDINARY DIFFERENTIAL
C EQUATIONS USING THE FOURTH ORDER STANDARD RUNGE-KUTTA METHOD.
C
C *PARAMETERS:
C
C DERIV- IS THE EXTERNAL NAME OF THE USER'S ROUTINE TO EVALUATE THE
C         RIGHT HAND SIDES OF THE DIFFERENTIAL EQUATIONS.  CALLED AS
C                 CALL DERIV(Y,YP,N,X)
C         WHERE:
C         Y   - IS A REAL ARRAY OF LENGTH "N" CONTAINING THE DEPENDENT
C                 VARIABLES Y(I).
C         YP  - IS A REAL ARRAY IN WHICH THE DERIVATIVES OF THE Y(I)
C                ARE RETURNED BY "DERIV"
C         N   - IS THE NUMBER OF EQUATIONS BEING INTEGRATED.
C         X   - IS THE INDEPENDENT VARIABLE.
C Y0   - IS THE INITIAL VALUE FOR Y(I). A REAL ARRAY OF LENGTH N
C YP0 -  IS THE INITIAL VALUE OF THE DERIVATIVES OF Y(I) AT X,
C         INPUT IF IX.NE.0 AND OUTPUT OTHERWISE.
C Y1   - IS THE ARRAY OF FINAL Y(I) VALUES AT X+H
C N    - IS THE NUMBER OF EQUATIONS BEING INTEGRATED.
C X    - IS THE INDEPENDENT VARIABLE RETURNED AS X+H
C H    - IS THE FIXED SINGLE STEP TAKEN IN X.
C IX   - IS ZERO IF THE DERIVATIVES AT THE INPUT VALUE X HAVE TO
C         BE CALCULATED, ELSE IS NON-ZERO AND THE DERIVATIVES ARE
C         INPUT IN YP0.
C WK   - IS A WORKING ARRAY OF AT LEAST 3*N REAL ELEMENTS.
C
C INPUT VARIABLES:  DERIV, Y0, N, X, H, IX AND POSSIBLY YP0
C OUTPUT VARIABLES: YP0, Y1, X
C
C NOTE.  THE VECTORS Y, YP0 AND Y1 MUST ALL BE DISTINCT - AND DISTINCT
C FROM THE COLUMNS OF WK!!
C
C *METHOD.
C 	THE FOURTH ORDER RUNGE-KUTTA ALGORITHM IS CARRIED IN THE
C STRAIGHTFORWARD WAY.  THE CODING IS DESIGNED TO REDUCE TO A MINIMUM
C THE NUMBER OF FLOATING POINT MULTIPLICATIONS.
C 	IF IX IS ZERO THE DERIVATIVES AT THE INITIAL X VALUE ARE CALCULATED
C IN YP0.  OTHERWISE THEY ARE ASSUMED TO BE THERE ALREADY.
C
C *ACCURACY.
C 	ROUNDING ERROR WILL BE OF THE ORDER OF THE MACHINE SIGNIFICANCE
C TRUNCATION ERRORS DEPEND ON THE STEP SIZE AND THE PROBLEM.
C
C *RESTRICTIONS.
C
C *ERROR CONDITIONS.
C
C *NON STANDARD ROUTINES CALLED.
C 	DERIV(USER'S)
C
C *TYPICAL TIMES.
C
C *ORIGIN.   M.R.MANNING.
C
C *COMMENTS.
C
C #END.
C
C
C #START.
      HBY2=H/2.
      HBY6=H/6.
C
      IF(IX) 10,5,10
5     CALL DERIV(Y0,YP0,N,X)
C
10    DO 12 I=1,N
12    Y1(I)= Y0(I) + HBY2*YP0(I)
      XA=X+HBY2
      CALL DERIV(Y1,WK,N,XA)
      DO 14 I=1,N
14    Y1(I)= Y0(I) + HBY2*WK(I,1)
      CALL DERIV(Y1,WK(1,2), N,XA)
      DO 16 I=1,N
16    Y1(I)= Y0(I) + H*WK(I,2)
      X=X+H
      CALL DERIV(Y1,WK(1,3), N,X)
      DO 20 I=1,N
      T=WK(I,1) + WK(I,2)
      T=HBY6*(YP0(I) + T + T +WK(I,3) )
      WK(I,1)=T
20    Y1(I)=Y0(I) + T
C
      RETURN
      END
����������������������������������������������������������������������������������������������������������������������������������������������������������������������������������������