Next: Discussion Up: The Fourth-Order Runge-Kutta Method Previous: Adaptive Stepsize Control

## Example Code

```SUBROUTINE rkqs(y, dydx, x, htry, eps, yscal, hdid, hnext, derivs)

USE nrtype; USE nrutil, ONLY : assert_eq, nrerror
USE nr, ONLY : rkck

IMPLICIT NONE

REAL(sp), DIMENSION(:), INTENT(inout) :: y
REAL(sp), DIMENSION(:), INTENT(in) :: dydx, yscal
REAL(sp), INTENT(inout) :: x
REAL(sp), INTENT(in) :: htry, eps
REAL(sp), INTENT(out) :: hdid, hnext

INTERFACE
SUBROUTINE derivs(x, y, dydx)
USE nrtype
IMPLICIT NONE
REAL(sp), INTENT(in) :: x
REAL(sp), DIMENSION(:), INTENT(in) :: y
REAL(sp), DIMENSION(:), INTENT(out) :: dydx
END SUBROUTINE derivs
END INTERFACE

INTEGER(i4b) :: ndum
REAL(sp) :: errmax, h, htemp, xnew
REAL(sp), DIMENSION(SIZE(y)) :: yerr, ytemp
REAL(sp), PARAMETER :: safety=0.9_sp, pgrow=-0.2_sp, pshrnk=-0.25_sp, &
errcon=1.89e-4 ! (5/safety)**(1/pgrow)

ndum=assert_eq(SIZE(y), SIZE(dydx), SIZE(yscal), 'rkqs')
h=htry
DO
CALL rkck(y, dydx, x, h, ytemp, yerr, derivs)
errmax=MAXVAL(ABS(yerr(:)/yscal(:)))/eps
IF (errmax <= 1.0) EXIT
htemp=safety*h*(errmax**pshrnk)
h = SIGN(MAX(ABS(htemp), 0.1_sp*ABS(h)), h)
xnew = x+h
IF (xnew == x) CALL nrerror('stepsize underflow in rkqs')
END DO
IF (errmax > errcon) THEN
hnext=safety*h*(errmax**pgrow)
ELSE
hnext=5.0_sp*h
END IF
hdid=h
x=x+h
y(:)=ytemp(:)

END SUBROUTINE rkqs
```

Zdzislaw Meglicki
2001-02-26