Let's consider differential equation: $$x^{(n)}(t) = f(t,x(t),x'(t),\ldots,x^{(n-1)}(t))$$ with initial conditions $x(0) = x_0, x'(0) = x_1, \ldots, x^{(n-1)}(0) = x_{n-1}$. It is easy to solve this equation numerically by substitution $X_1 = x, X_2 = x',\ldots, X_{n-1} = x^{(n-1)}$ and using Euler or Runge-Kutta recurence schemes.
The analogous method cannot be applied to the equation: $$D^{\alpha} x(t) = g(t,x(t),D^{\beta}x(t))$$ where $D^{\alpha}, D^{\beta}$ are the Caputo derivatives, $0<\beta<\alpha$, and initial conditions are $x(0) = x_0, x'(0) = x_1, \ldots, x^{(n-1)}(0) = x_{n-1}$, where $n = \lceil \alpha \rceil$. How to solve such equations numerically? (if $g$ depends only on $t$ and $x(t)$, then I can solve this equation,for example by using Grunwald-Letnikow derivative, the problem is second derivative with non-integer order and nonlinear terms like $x(t) \cdot D^{\beta} x(t)$).