This algorithm is ported for Java from C++ code (last updated: Nov 15, 2002), written by Blake Ashby on the base of code RADAU5 originally written in FORTRAN (version of July 9, 1996, latest small correction: January 18, 2002) by E. Hairer and G. Wanner (Universite de Geneve, Dept. de Mathematiques Ch-1211 Geneve 24, Switzerland) and which is in details described in their book[1].
This code computes the numerical solution of a differential system of first order ordinary differential equations `y'=f(x,y).`

The method used is an implicit Runge-Kutta method (RADAU IIA) of order 5 with step size control and continuous output. (See Section IV.8 of [1]). The method is designed to solve stiff systems, for non-stiff systems it is rationally to use simpler methods (requiring less computational load), e.g. Dormand-Prince.

User-provided solver parameters Description

itoler

itoler = 0: Both rtoler and atoler are scalars. The code keeps, roughly, the local error of y[i] below ` rtoler*fabs(y[i])+atoler `.
itoler = 1: Both rtoler and atoler are vectors. The code keeps the local error of y[i] below ` rtoler[i]*fabs(y[i])+atoler[i] `.

hinit

Initial step size guess;
For stiff equations with initial transient, ` h = 1.0/(norm of f') `, usually 1.0e-3 or 1.0e-5, is good. This choice is not very important, the step size is quickly adapted. (if ` h = 0.0 ` on input, the code sets ` h = 1.0e-6 `).

hmax

Maximal step size, default (when hmax is set = 0) ` xend - x `.

nmax

This is the maximal number of allowed steps. The default value is 100000.

uround

rounding unit, default 1.0e-16.

safe

The safety factor in step size prediction, default 0.9.

facl
facr

Parameters for step size selection the new step size is chosen subject to the restriction ` 1/facl <= hnew/hold <= 1/facr `.
Default values: ` facl = 5.0 `, ` facr = 1.0/8.0 `.

nit

The maximum number of Newton iterations for the solution of the implicit system in each step. The default value is 7.

startn

If ` startn == 0 ` the extrapolated collocation solution is taken as starting value for Newton's method. If ` startn != 0 ` starting values are used.
The latter is recommended if Newton's method has difficulties with convergence. (This is the case when nstep is larger than ` naccpt + nrejct `; see output parameters). Default is ` startn = false `.

npred

Switch for step size strategy:

If ` npred == 1 `: mod. predictive controller (Gustafsson)
If ` npred == 2 `: classical step size control.

The default value (for ` npred = 0 `) is ` npred = 1 `.
The choice ` npred = 1 ` seems to produce safer results. For simple problems, the choice ` npred == 2 ` produces often slightly faster runs.

hess

If ` hess != 0 `, the code transforms the Jacobian matrix to Hessenberg form. This is particularly advantageous for large systems with full Jacobian. It does not work for banded Jacobian (` mljac < n `) and not for implicit systems (` imas = 1 `).

fnewt

Stopping criterion for Newton's method, usually chosen < 1. Smaller values of fnewt make the code slower, but safer.
Default ` min(0.03, sqrt(rtoler)) `.

quot1
quot2

If ` quot1 < hnew/hold < quot2 `, then the step size is not changed. This saves, together with a large thet, lu-decompositions and computing time for large systems.
For small systems one may have ` quot1 = 1.0 `, ` quot2 = 1.2 `, for large full systems ` quot1 = 0.99 `, ` quot2 = 2.0 ` might be good. Defaults ` quot1 = 1.0 `, ` quot2 = 1.2 `.

thet

Decides whether the Jacobian should be recomputed. Increase thet, to 0.1 say, when Jacobian evaluations are costly. for small systems thet should be smaller (0.001, say). Negative thet forces the code to compute the Jacobian after every accepted step. Default 0.001.

## References

1. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems. Springer Series in Computational Mathematics 14, Springer-Verlag 1991, Second Edition 1996.