万有引力の運動方程式を解くと楕円軌道になることはわかっているのだが，楕円軌道は $r(\phi)$ で表され，時間 $t$ の陽関数としては表せていないので，$\phi$ を $t$ の関数として数値的に求めるための下ごしらえ。
「万有引力の運動方程式を数値的に解く準備」では，万有引力の運動方程式を数値的に解く前の下ごしらえをしていた。実際に Python で解くのは「Python で万有引力の運動方程式を数値的に解いて楕円軌道のアニメをつくる」に書いている。
連立2階常微分方程式としての（規格化された）運動方程式を数値的に解くのも教育的ではあるが，この問題は楕円軌道になることがわかっている。ただし，解が
$$r = \frac{a (1 – e^2)}{1 + e \cos\phi}$$
と書けることはわかっているのだが，$r$ も $\phi$ も時間 $t$ の関数としてはあからさまに与えられていないため，ある時刻に楕円軌道上のどこにいるかは，数値的に解いてみないとわからない。
ここでは，楕円軌道は $r(\phi)$ で表されるので，$\phi$ を $t$ の関数として数値的にもとめ，$\phi(t)$ が数値的に求まれば，
$$x(t) = r(\phi(t))\,\cos\phi(t), \quad y(t) = r(\phi(t))\, \sin\phi(t)$$
のように，時刻を与えれば楕円軌道上の位置がわかるようにしてみる。以下，「参考：ニュートン力学における万有引力の2体問題」から要点を抜粋して利用する。
万有引力の2体問題は1体問題に帰着する
万有引力の2体問題は1体問題に帰着し，運動方程式は
$$\frac{d^2\boldsymbol{r}}{dt^2} = – \frac{GM \boldsymbol{r}}{r^3}$$
保存量（運動の定数）
運動方程式から，2つの保存量（運動の定数）が得られる。1つは，(単位質量あたりの）角運動量の $z$ 成分に相当する $\ell$：
$$\ell \equiv r^2 \frac{d\phi}{dt}$$
もう1つは（単位質量あたりの）力学的エネルギーに相当する $\varepsilon$：
$$\varepsilon \equiv \frac{1}{2} \left\{\left(\frac{dr}{dt}\right)^2 + \frac{\ell^2}{r^2}\right\} – \frac{GM}{r}$$
解：楕円軌道
$$r = \frac{a(1-e^2)}{1 + e \cos\phi}$$
ここで，軌道長半径 $a$，離心率 $e$ と保存量の関係は
$$\ell = \sqrt{GM a (1-e^2)}, \quad \varepsilon = – \frac{GM}{2a}$$
解くべき方程式
基本的に以下の1階常微分方程式を解くだけでよい。
$$\frac{d\phi}{dt} = \frac{\ell}{r^2}, \quad r = \frac{a(1-e^2)}{1 + e \cos\phi}$$
無次元化
ケプラーの第3法則
$$\frac{a^3}{P^2} = \frac{GM}{4\pi^2}, \quad \therefore\ \ P = \frac{2\pi a^2}{\sqrt{GM a}}$$
軌道周期 $P$ を使って無次元化された時間 $T$ を
$$T \equiv \frac{t}{P}$$
とすると，
$$\frac{d\phi}{dT} = 2 \pi \frac{(1+e \cos\phi)^2}{(1-e^2)^{3/2}} \tag{1}$$
これを数値的に解き，$\phi(T)$ がわかったら，
\begin{eqnarray}
R(T) &\equiv& \frac{r}{a} = \frac{1-e^2}{1+e\cos\phi(T)} \tag{2}\\
X(T) &\equiv& R(T) \cos\phi(T) \tag{3}\\
Y(T) &\equiv& R(T) \sin\phi(T) \tag{4}
\end{eqnarray}
として，規格化された $X$ 座標，$Y$ 座標を求めればよい。