This is a supplemental answer for now because while we know that a two body orbit can be reduced to a one body orbit around a central potential, doing that here will be a little distracting and I think the result for the one body in central potential looks cleaner. See also answers to Can the radial oscillations of an elliptical orbit be solved using a fictitious centrifugal potential?
Per this comment I know I've had a discussion somewhere in this site (or in Astronomy SE) where it was first explained to me that Kepler orbits do have analytical solutions you can write down for time as a function of position, even though we still do need to use numerical techniques (e.g. Newton's method) to solve position as a function of time. (see also How did Newton and Kepler (actually) do it?)
If someone finds it before I do please feel free to add a link here, thanks!
Equation 27 in Wikipedia's Kepler orbit; Properties of trajectory equation is
$$t = a \sqrt{\frac{a}{\mu}}\left(E - e \sin E \right)$$
where $a$ is the semimajor axis, $\mu$ is the standard gravitational parameter also known as the product $GM$, $e$ is the eccentricity and $E$ is the Eccentric anomaly.
The relationship between $E$ and the true anomaly $\theta = \arctan2(y, x)$ is
$$\tan \frac{\theta}{2} = \sqrt{ \frac{1+e}{1-e} } \tan \frac{E}{2}$$
and solving for $E$:
$$E(\theta) = 2 \arctan \sqrt{ \frac{1-e}{1+e} } \tan \frac{\theta}{2}.$$
plugging back in to the first equation (but not writing it all out):
$$t(\theta) = a \sqrt{\frac{a}{\mu}}\left(E(\theta) - e \sin E(\theta) \right)$$
Let's try a numerical check of this amazing result. Note that with $a=1$ and $\mu=1$ the period is $2 \pi$.
The last plot at bottom left shows that the analytical $t(\theta)$ based on $\theta$ from a numerically integrated orbit matches the time used in the numerical calculation for an $e=0.8$ elliptical orbit. There will be numerical glitches or singularities at the endpoints and for $e=1$ but it seems to check out nicely!
Python script:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
def deriv(X, t):
x, v = X.reshape(2, -1)
acc = -x * ((x**2).sum())**-1.5
return np.hstack((v, acc))
halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
e = 0.8
a = 1.0
mu = 1.0
r_peri, r_apo = a*(1.-e), a*(1.+e)
v_peri, v_apo = [np.sqrt(2./r - 1./a) for r in (r_peri, r_apo)]
T = twopi * np.sqrt(a**3/mu)
X0 = np.array([r_peri, 0, 0, v_peri])
X0 = np.array([-r_apo, 0, 0, -v_apo])
times = np.linspace(-T/2., T/2., 1001)
answer, info = ODEint(deriv, X0, times, full_output=True)
x, y = answer[1:-1].T[:2]
theta = np.arctan2(y, x)
E = 2. * np.arctan(np.sqrt((1.-e)/(1.+e)) * np.tan(theta/2))
t = a * np.sqrt(a/mu) * (E - e * np.sin(E))
if True:
plt.figure()
plt.subplot(2, 1, 1)
plt.plot(x, y)
plt.plot([0], [0], 'ok')
plt.gca().set_aspect('equal')
plt.title('y vs. x numerical')
plt.subplot(2, 1, 2)
plt.plot(times[1:-1], x)
plt.plot(times[1:-1], y)
plt.xlim(-pi, pi)
plt.title('x(t) and y(t) numerical')
plt.show()
plt.subplot(2, 2, 1)
plt.title('theta(t_numerical)')
plt.plot(times[1:-1], theta)
plt.xlim(-pi, pi)
plt.ylim(-pi, pi)
plt.gca().set_aspect('equal')
plt.subplot(2, 2, 2)
plt.title('E_analytic(theta_numerical)')
plt.plot(E, theta)
plt.xlim(-pi, pi)
plt.ylim(-pi, pi)
plt.gca().set_aspect('equal')
plt.subplot(2, 2, 3)
plt.title('theta(t_analytic)')
plt.plot(t, theta)
plt.xlim(-pi, pi)
plt.ylim(-pi, pi)
plt.gca().set_aspect('equal')
plt.subplot(2, 2, 4)
plt.title('t_analytic(t_numerical)')
plt.plot(t, times[1:-1])
plt.xlim(-pi, pi)
plt.ylim(-pi, pi)
plt.gca().set_aspect('equal')
plt.show()
