Modeling the distance between two point masses in space with gravity as a function of time

(NOTE: When I started writing this response I had no idea it was going to become so long. At some point, I decided I had put in too much time to not finish it.)

I don't think this problem has an analytic solution for $R(t)$, but perhaps I'm wrong. Your expression for $ddot{R}$ certainly looks correct, so I'm not sure where the problem is. Here is another approach.

Kepler's third law will still apply to this system:

$$tau^2 = frac{4pi^2a^3}{Gleft(m_1 + m_2right)}.$$

For infalling particles, the period $tau$ is twice the time $t_{textrm{col}}$ needed for collision, and the semi-major axis $a$ is twice the starting separation $r_0$. So we may write

$$t_{textrm{col}} = pisqrt{frac{r_0^3}{8GM}},$$

where $M = left(m_1 + m_2right)$.

In the center-of-mass frame, $m_1r_1 = m_2r_2$. If we define $rho equiv r_1 + r_2$, we can show after a few lines of algebra that

$$murho = m_1r_1 = m_2r_2,$$

where $mu$ is the reduced mass:

$$mu equiv frac{m_1m_2}{m_1 + m_2}.$$

Differentiate each term in the $murho$ expression to get

$$mudot{rho} = m_1dot{r_1} = m_2dot{r_2}.$$

We can identify $dot{rho}$ as $v$, $mu v$ as the momentum of each particle, and $frac{1}{2}mu v^2$ as the total kinetic energy. Since mechanical energy is conserved, we may write

$$frac{1}{2}mu v^2 - frac{Gm_1m_2}{rho} = 0 - frac{Gm_1m_2}{rho_0}.$$

This leads to

$$v^2 = frac{2Gm_1m_2}{mu}left(frac{1}{rho} - frac{1}{rho_0}right).$$

Simplifying,

$$v = sqrt{2Gleft(m_1 + m_2right)left(frac{1}{rho} - frac{1}{rho_0}right)} = sqrt{2GMleft(frac{1}{rho} - frac{1}{rho_0}right)}.$$

The usual approach from here is to introduce dimensionless quantities:

$$eta equiv frac{rho}{rho_0}, textrm{ and}$$ $$psi equiv frac{t}{t_{textrm{col}}},$$

where $t_{textrm{col}}$ is now written as $pisqrt{rho_0^3/8GM}$. The substitutions become

$$rho = rho_0eta implies drho = rho_0 deta, textrm{ and}$$

$$t = t_{textrm{col}}psi implies dt = t_{textrm{col}}dpsi.$$

Since $v = -drho/dt$, we get

$$v = -frac{rho_0}{t_{textrm{col}}}frac{deta}{dpsi}.$$

Substituting for $t_{textrm{col}}$,

$$v = -frac{rho_0}{pisqrt{rho_0^3/8GM}}frac{deta}{dpsi} = -frac{1}{pisqrt{rho_0/8GM}}frac{deta}{dpsi}.$$

Substituting into the expression for $v$ above,

$$-frac{1}{pisqrt{rho_0/8GM}}frac{deta}{dpsi} = sqrt{2GMleft(frac{1}{etarho_0} - frac{1}{rho_0}right)},$$

$$frac{deta}{dpsi} = -pisqrt{frac{rho_0}{8GM}frac{2GM}{rho_0}left(frac{1}{eta} - 1right)},$$

$$frac{deta}{dpsi} = -frac{pi}{2}sqrt{left(frac{1}{eta} - 1right)}.$$

Separate the variables to get

$$frac{pi}{2}dpsi = frac{-deta}{sqrt{left(frac{1}{eta} - 1right)}}.$$

The next step is to integrate both sides. In the center-of-mass frame coordinates, the particles fall from an initial separation distance $rho_0 = rho$ at $t = 0$ to $rho = 0$ at $t = t_textrm{col}$. So $eta$ varies from an initial value of 1 to a final value of 0, while $psi$ varies from an initial value of 0 to a final value of 1.

The integrals for intermediate values of $eta$ and $psi$ are therefore

$$frac{pi}{2}int_{0}^{psi}dpsi = int_{1}^{eta}frac{-deta}{sqrt{left(frac{1}{eta} - 1right)}}.$$

At this point, we can switch notation to something that is more familiar. Trading $psi$ for $t$ and $eta$ for $r$ (just be very careful going backward from here to check the steps!), we can write

$$tfrac{pi}{2} = int_{1}^{r}frac{-dr}{sqrt{left(frac{1}{r} - 1right)}}.$$

There are probably several ways to solve the $r$ integral. I used an online tool:

$$tfrac{pi}{2} = rsqrt{frac{1}{r}-1} + tan^{-1}left(sqrt{frac{1}{r} - 1}right).$$

The $tan^{-1}$ function can probably be simplified (I have seen $cos^{-1}$ in several sources), but I don't think there is any way to invert this to get an expression for $r$ as a function of $t$. The best we could do is get numerical solutions.