楕円軌道は $r(\phi)$ で表されるので，$\phi$ を $t$ の関数として数値的に求めるための下ごしらえ。
「万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ」では，万有引力の2体問題（は結局1体問題に帰着するんだけど）の運動方程式を数値的に解く前の下ごしらえをしていた。実際に Maxima で解くのは「下ごしらえした万有引力の2体問題の運動方程式を Maxima で数値的に解く」に書いている。
連立2階常微分方程式としての（規格化された）運動方程式を数値的に解くのも教育的ではあるが，この問題は楕円軌道になることがわかっている。ただし，解が
$$r = \frac{a (1 – e^2)}{1 + e \cos\phi}$$
と書けることはわかっているのだが，$r$ も $\phi$ も時間 $t$ の関数としてはあからさまに与えられていないため，ある時刻に楕円軌道上のどこにいるかは，数値的に解いてみないとわからない。
時刻 $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$ 座標を求めればよい。