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