# Algorithm

On the first step, CVODE(S) will use backward Euler, since both Adams-Moulton and BDF reduce to this at order 1. For this step, the local truncation error (LTE) has the form

(1)We search for an initial step $h$ that roughly solves

(2)For this, we must estimate $\ddot y_0 = \ddot y (t_0)$, knowing $\dot y_0 = f(t_0, y_0)$. For a given stepsize $\bar h$, we use the approximation

(3)and iterate until the condition (2) is satisfied.

To get started, we use bounds for $|h|$ based on the initial time $t_0$ and the first requested output time $t_{out}$. The lower bound is set to

(4)where $\epsilon$ is the unit roundoff, while the upper bound is set to

(5)possibly adjusted downward to ensure that

(6)We start the iterations with $\bar h = \sqrt{h_L h_U} \cdot \text{sign} (t_{out} - t_0 )$ and we stop when $1/2 < h/\bar h < 2$ (this is good enough, since the error in the first step may reset the initial stepsize).

In CVODES, the above algorithm is extended to include the quadrature (if quadrature variables exist and they are included in the error test) and sensitivity variables (if forward sensitivity analysis is enabled and the sensitivity variables are included in the error test). In such a case, all norm evaluations (the max-norm in estimating $h_U$ and the wrms-norm of the second derivative) are replaced with the maximum of the corresponding norms over all states, quadrature variables, and sensitivity variables. For example, under the appropriate conditions, the stopping criteria (2) may be replaced with

(7)where $y^Q$ is the vector of quadrature variables and $y^{S_i}$ are the $N_s$ sensitivity vectors.

# Implementation details

In attempting to find an acceptable initial stepsize, CVODE(S) will iterate a maximum MAX_ITERS (=4) times after which it will return (the resulting stepsize may be later reset depending on the error in the first step). Moreover, the iterations will be stoped and a *fall back* value ($\bar h$) will be used if $hnew/hg > 2$ after one iteration as this probably means that the $\ddot y_0$ value is bad because of cancellation error.

The computation of the initial stepsize will return a failure flag if

- $t_{out}$ is too close to $t_0$; i.e. $|t_{out} - t_0| < 2 \cdot \epsilon \cdot \max \{ |t_0| , |t_{out}| \}$ (CVode return flag: CV_TOO_CLOSE)
- a user-supplied function (ODE RHS or, if applicable, quadrature or sensitivity RHS) failed unrecoverably; i.e. it returned with a negative value (CVode return flag: one of CV_*RHSFUNC_FAIL)
- a user-supplied function failed recoverably, i.e. it returned a positive value, more than MAX_ITERS times (CVode return value: one of CV_REPTD_*RHSFUNC)