where \(\mathbf {y}^\text{BDF}_j = \mathbf {y}^\text{BDF}(t_{j})\) . The order \(q\)  can vary between 1 and 5, and is adapted simultaneously with the step size \(h_j\) . The coefficients \(\alpha _{j,i}\) , \(\beta _{j}\)  are determined by the recent history of the step sizes [1]}, [2]}. For this method, estimation of the local truncation error and adaptation rules are highly complex and all details are described in [3]}. A nonlinear rootfinding problem must be solved at each step to obtain \(\mathbf {y}^\text{BDF}_j\) , and this is done approximately using modified Newton iteration. This in turn requires solving a linear system, for which Stan uses the dense linear solver facilities of SUNDIALS. In addition to controlling the estimated local truncation error, also the Newton iteration convergence has to be controlled. The stopping criterion for the Newton iterations is in CVODES attempted to set so that its error does not interfere with local error control.
