Discretization of ODE system I am fairly new to the discretization of ODE systems (indeed a good reference would be helpful). I have a system of ODEs that basically looks like this
$$
\begin{align}
    \frac{d x(t)}{dt} &= v(t) \\
    \frac{d v(t)}{dt} &= a(t, x_t, v_t)
\end{align}
$$
How do I discretize this and , given a discretization, how do I know if it's correct? I have seen in an example online something like this:
$$
\begin{align}
    x_{t + \delta/2} &= x_t + \frac{\delta}{2} v_t \\
    v_{t + \delta/2} &= v_t + \frac{\delta}{2}a(t, x_{t + \delta/2}, v_t)\\
    x_{t + \delta }  &= x_{t + \delta/2} + \frac{\delta}{2} v_{t + \delta/2}
\end{align}
$$
Does it make sense? is this a correct discretization and why? what properties does it have?
 A: There is indeed a strategy in the approximation of second order DE or PDE like the wave equation that uses a dual or leapfrogged lattice. This is used to implement approximations of the derivatives stack as central difference quotients without spreading the support of the formulas too far
$$
\frac{y_{t+δ}-y_t}δ\approx y'_{t+δ/2},\\
\frac{y'_{t+δ/2}-y'_{t-δ/2}}δ\approx y''_t.
$$
If the formula for the right side does not satisfy the graduation of the leapfrogged scheme, as here for $y''=a(t,y,y')$, one can repair that like in the midpoint method
$$
y''_t\approx a\left(t,y_t,\frac{y'_{t+δ/2}+y'_{t-δ/2}}2\right)
$$
or as in the trapezoidal method
$$
y''_t\approx\frac{a(t,y_t,y'_{t+δ/2})+a(t,y_t,y'_{t-δ/2})}2
$$
The next challenge is that these will then turn out to give implicit formulas for $y'_{t+δ/2})$. These get usually solved iteratively. One can make the iteration count fixed to get an explicit method à la the explicit Heun method, or with a variable number of steps until the error in the implicit formula is sufficiently reduced to be lower than the truncation error of the method.

In general this scheme will give no further advantages over its being of second error order and its compact support. In the case where $My''=-\nabla V(y)$ is the equation of motion for the Hamiltonian/energy function $H(y,y')=\frac12y'^TMy'+V(y)$, that is, $a(t,y,y')=-M^{-1}∇V(y)$, this leapfrogged method gives a variant of the Verlet method. This is still a second order method in general, but the energy oscillates inside a band of width $O(δ^2)$, which stays invariant with order 4 or better. The phase of the solution however will usually have a second order error.
