This code is really just a wrapper.

The variables are similar to the ones used in the previous Runge Kutta subroutine:

**y**- the equation that is being solved is of the form

so this variable represents a vector of*y*_{k}values, where*k*runs through the equations. On entry**y**stands for . On exit it will stand for . **dydx**- is the initial value of
derivs(
*x*,**y**), that is . **x**- is
*x*_{n} **htry**- is the initial guess for a good that may be changed if a requested accuracy of integration is not met.
**eps**- The error associated with the whole system
of equations is going to be evaluated in the following
manner. The Runge-Kutta subroutine
`rkck`

to be called by`rkqs`

will return a new value plus an absolute error estimate . That error is then going to be*scaled*by dividing by a user supplied array`yscal`

. This new*scaled*error is then going to be compared to , which is what this parameter`eps`

stands for. **yscal**- is the scaling array. If you're happy with
unscaled values of
**y**simply set it to 1. **hdid**- is the value of that has really been used, after all the mucking up, to make the step.
**hnext**- is the suggested next value of .
**derivs**- is the subroutine used to evaluate the right hand
side of

As most other codes discussed in this lecture notes
the action here begins by checking that the dimensions
of variables passed to the subroutine are correct. The
dimension of
**y** that is extracted is
discarded:

ndum=assert_eq(SIZE(y), SIZE(dydx), SIZE(yscal), 'rkqs')

Then the step size
is set to the suggested
step size of `htry`

:

h=htryand we enter the

`DO`

loop within which we
- 1.
- call
`rkck`

to evaluate the next**y**for the suggested value of , but also to give us the truncation error in`yerr`

:CALL rkck(y, dydx, x, h, ytemp, yerr, derivs)

- 2.
- scale the returned error and compare it to :
errmax=MAXVAL(ABS(yerr(:)/yscal(:)))/eps IF (errmax <= 1.0) EXIT

observe that since we are working with the system of equations here, and each equation is going to have its own different value of error , we pick up the largest error and make that stand for the error for the whole system of equations. If the error is within the prescribed limits we accept the result and exit the`DO`

loop. - 3.
- If the error is too large, then we have to repeat
the step with a shorter .
So we shrink the
step, but by no more than a factor of 10:
htemp=safety*h*(errmax**pshrnk) h = SIGN(MAX(ABS(htemp), 0.1_sp*ABS(h)), h)

- 4.
- then check if the step size hasn't shrunk too much, in
which case the subroutine aborts with an error message:
xnew = x+h IF (xnew == x) CALL nrerror('stepsize underflow in rkqs')

- 5.
- and return to the top of the loop when the next trial Runge Kutta step is made.

After we have finished with the looping and have the new
values of
**y** as well as the new value of error,
we suggest stretching ,
but by no more than 5 times:

IF (errmax > errcon) THEN hnext=safety*h*(errmax**pgrow) ELSE hnext=5.0_sp*h END IF

The last three lines move the used
into `hdid`

,
the new value of
**y** goes into `y`

, and
*x* itself becomes updated to
:

hdid=h x=x+h y(:)=ytemp(:)

The subroutine `rkck`

is going to be very similar to
subroutine `rk4`

, but instead of making just one
Runge-Kutta step, it has to make

- 1.
- one full size step
- 2.
- two half size steps

`ytemp`

and the observed error
in `yerr`

.
I leave it to you to develop your own version of `rkck`

.
The easiest way to do that is to rewrite subroutine `rk4`

,
which we have discussed in section 4.4.1.