I will phrase the question only for second order ODE's although the same question might come up for higher orders.
ODEs of the form
$$ f(t, x,\dot{x}) = \ddot{x} $$
are usually converted to first order ODEs with the trick $y=(x,\dot{x})$ which then results in an equivalent formulation
$$ g(t, y) = \begin{pmatrix} y_2 \\ f(t, y) \end{pmatrix} =\dot{y} $$
But while discretizations of first order ODE's are well studied with arbitrary convergence order being achieved by the family of Runge Kutta methods, I am not aware of any theory on higher order ODE discretizations. And I believe that the conversion to first order ODE's probably sacrifices a lot of potential (I will explain this belief in my motivation).
Motivation
We will be looking at gradient decent with momentum on a loss function $L$. Instead of setting the velocity of the decent equal to the steepest decent
$$ \dot{x} = -\nabla L(x) $$
we instead want to set the acceleration to be equal to the steepest decent and add some decelerating force ("friction") proportional to the current velocity
$$ \ddot{x} = -\nabla L(x) - \alpha\dot{x} $$
The trick above and euler discretization yields
$$ y^{n+1} = y^n + \eta \begin{pmatrix} y^n_2 \\ -\nabla L(y^n_1) - \alpha y^n_2 \end{pmatrix} $$
with a bit of reparametrization ($h=\eta^2$, $x=y_1$, $m=y_2/\eta$, $\beta=1-\eta\alpha$) this (almost) results in the well known heavy ball (momentum) method
$$ \begin{aligned} x^{(n)} &= x^{(n-1)} + hm^{(\color{red}{n-1})}\\ m^{(n)} &= \beta m^{(n-1)} - \nabla L(x^{(n-1)}) \end{aligned} $$
In reality (e.g. https://distill.pub/2017/momentum/) people use the current momentum $m^{(n)}$ for the position update. Now taking the learning rate to zero ($h\to 0$), this will likely result in the same ODE. But in the discrete case it makes a difference.
Substituting in $m^{(n-1)}$ into the first equation we would get
$$ x^{(n)} = x^{(n-1)} + h(\beta m^{(n-1)} - \nabla L(x^{(n-2)}) $$ which means we are not using the gradient information from our previous position $x^{(n-1)}$ to calculate the current position $x^{(n)}$ but the one before that. This means our optimization algorithm will overshoot when a change in gradient occurs.
Nesterov's momentum goes one step further and first applies the current momentum and only then calculates the gradient
$$ x^{(n)} = \underbrace{x^{(n-1)} + h\beta m^{(n)}}_{\text{"momentum move"}} -\nabla L(x^{(n-1)} + h\beta m^{(n)}) $$
This is argued to perform even better (cf. https://stats.stackexchange.com/a/368179/265966)
Whatever jump comes first, my Momentum Jump would be the same. So I should consider the situation as if I have already made my Momentum Jump, and I am about to make my Slope Jump.
Concluding Thoughts
Now all these slight heuristic modifications make intuitive sense. But I wonder if there is some theory to put this on a solid foundation and maybe even find an "optimal" discretization.
What I find remarkable is, that higher order methods (like the Heun's method) also do some form of look ahead:
$$ x^{(k+1)} = x^{(k)} + \tfrac12 h \big[f(x^{(k)}) + f\big(\underbrace{x^{(k)} + hf(x^{(k)})}_{\text{"look ahead"}}\big)\big] $$
In our second order ODE we just did a skewed lookahead, where we used a newer momentum to update our position. So higher order ODE discretization might have to somehow incorporate that.
