$$\begin{aligned}
\beta &=& \frac{2M}{r_0}-v_0^2\\
\eta_0 &=& \mathbf{r}_0 \cdot \mathbf{v}_0 \\
\zeta_0 &=& M-\beta r_0\end{aligned}$$ where $r_0 = |\mathbf{r}_0|$ and $v_0 = |\mathbf{v}_0|$. Note that the semi-major axis $a$ can be written as $a=M/\beta$.

Second, we need to solve Kepler’s equation which, using the above notation, takes the form $$\begin{aligned}
r_0 X + \eta_0 G_2(\beta,X) + \zeta_0 G_3(\beta,X) - \mathit{dt} = 0. \label{eq:kepler} \end{aligned}$$ We solve this equation for $X$. This is a non-algebraic (i.e. transcendental) equation that we need to solve iteratively, for example using Newton’s method. In Sect. \[sec:newton\], we describe our algorithm in detail.

Third, having solved Kepler’s Equation, we can calculate the so called Gau[ß]{} $f$ and $g$-functions as well as their time derivatives via $$\begin{aligned}
f &= 1 - M \frac{G_2}{r_0}\quad & \dot f &= -\frac{M\,G_1}{r_0r} \label{eq:fg}&\\
g &= dt - M G_3            & \dot g &= 1-\frac{M\,G_2}{r},\label{eq:fg2}\end{aligned}$$ where $r = r_0 + \eta_0 G_1 + \zeta_0 G_2$. Note that all the $G$-functions depend on $\beta$ and the $X$ value found in the second step.

Fourth, we write the final positions and velocities as a linear transformation of the initial conditions using the Gau[ß]{} $f$ and $g$-functions: $$\begin{aligned}
\mathbf{r} &= f \mathbf{r}_0 + g \mathbf{v}_0
&\mathbf{v} 