I'm trying to model the motion of one mass $m$ about another larger mass $M$.
I have the following relation:
\begin{align} \ddot x_0 &= -GMx_0/r^3, \\ \ddot x_1 &= -GMx_1/r^3, \end{align} where $G$ is the gravitational constant and $r = (x_0^2 + x_1^2)^{1/2} = ||x||$ is the distance of the orbiting body from the center of mass of the two bodies.
Now, to solve computationally, I want to express the above relations as a system of first order ODEs, so that I can use computational methods (e.g. Matlab or Scipy's integrate.solve_ivp) to find a solution.
However I'm not sure how to set this up when we have two equations that depend on each other. I know the process for equations with a single variable is to re-write the system by defining something like $y_1(t) = x_0'$ and $y_2(t) = x_0'' = {-GMx_o/r^{3}}''$.
How can I put this into a system of first order ODEs that can be solved via Runge-Kutta?
This has been asked before.
The two-body problem can be solved almost completely in closed form. The equations of motion can be integrated to obtain the equations of the orbits, which are always conic sections (Kepler's first law generalised to include the unbounded cases of the parabola and hyperbola). However, to express the motions of the bodies as function of time, it's necessary to solve Kepler's equation, for which there is no known simple solution in closed form.