C <<<<<<<<<<<<<<<<<<<<<<<< DERIVS - DERIVS - DERIVS >>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINE DERIVS(X, Y, DYDX, N, SUBR, ERRMSG) C C Calculate the derivative of Y at X for N varaibles and store them in DYDX. C IMPLICIT REAL*8 (A-H, O-Z) DIMENSION Y(1), DYDX(1) CHARACTER ERRMSG*80, SUBR*20 C C Check input. C ERRMSG = ' ' IF (X .EQ. 0.D0) THEN ERRMSG = 'Attempted division by zero in DERIVS, called from ' 1 // SUBR RETURN ENDIF C DYDX(1) = -Y(2) DO 10 I=2,N DYDX(I) = Y(I-1) - (DBLE(I-1)/X)*Y(I) 10 CONTINUE C RETURN END C <<<<<<<<<<<<<<<<<<<<<<<<<<<< RK4 - RK4 - RK4 >>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINE RK4(Y, DYDX, N, X, H, YOUT, SUBR, ERRMSG) C C Given values for N variables Y and their derivatives DYDX known as X, use C the fourth-order Runge-Kutte method to advance the solution over an interval C H and return the incremented variables as YOUT, which need not be a distinct C array from Y. The user supplies the subroutine DERIVS(X, Y, DYDX) which C returns derivatives DYDX at X. C C Set the maximum number of functions. C PARAMETER (NMAX=10) IMPLICIT REAL*8 (A-H, O-Z) DIMENSION Y(N), DYDX(N), YOUT(N), YT(NMAX), DYT(NMAX), DYM(NMAX) CHARACTER ERRMSG*80, SUBR*20 C ERRMSG = ' ' C C Check number of functions that were passed. C IF (N .GT. NMAX) THEN ERRMSG = 'Number of functions passed exceeds NMAX in RK4, ' 1 'called from ' // SUBR RETURN ENDIF C HH = H*0.5D0 H6 = H/6.D0 XH = X+HH C C First step. C DO 11 I=1,N YT(I) = Y(I) + HH*DYDX(I) 11 CONTINUE C C Second step. C CALL DERIVS(XH, YT, DYM, 'RK1 - 2nd step', ERRMSG) IF (ERRMSG(0:0) .NE. ' ') RETURN DO 12 I=1,N YT(I) = Y(I) + HH*DYT(I) 12 CONTINUE C C Third step. C CALL DERIVS(XH, YT, DYM, 'RK1 - 3rd step', ERRMSG) IF (ERRMSG(0:0) .NE. ' ') RETURN DO 13 I=1,N YT(I) = Y(I) + H*DYM(I) DYM(I) = DYT(I) + DYM(I) 13 CONTINUE C C Fourth step. C CALL DERIVS(X+H, YT, DYT, 'RK1 - 4th step', ERRMSG) IF (ERRMSG(0:0) .NE. ' ') RETURN DO 14 I=1,N YOUT(I) = Y(I) + H6*(DYDX(I) + DYT(I) + 2.D0*DYM(I)) 14 CONTINUE C RETURN END