EDIT: I ran the $y'(t)$ through Wolfram|Alpha with $p = 1$, $v = 1$ and $k = -1$. The missing answer now comes out in terms of the arctangent function. So what is going on here? Why does it come out in terms of the arctangent instead of natural logarithm?
For $y(t)$ where $y$ is the distance between two point masses, classical mechanics gives us $y''(t) = \frac{k}{y(t)^2}$ for gravitational or electromagnetic attraction.
Multiplying each side by $y'(t)$ and integrating: $\int y''(t) y'(t) = \int \frac{k y'(t)}{y(t)^2}$ $\frac1{2}y'(t)^2 = c_1 - \frac{k}{y(t)}$
Initial position $p = y(0)$ and velocity $v = y'(0)$: $c_1 = \frac1{2}v^2 + \frac{k}{p}$
This brings us to: $y'(t) = \pm \sqrt{ v^2 + 2k\frac1{p} - 2k\frac1{y(t)}}$
Solving this differential equation results in an implicit function, basically the inverse of $y(t)$: $c_2 \pm t = \dfrac{k p^{3/2}}{(2k + pv^2)^{3/2}} \ln\left(y(t) \left(\sqrt{(2k - 2k\frac{p}{y(t)} + pv^2)(2k + pv^2)} + 2k + pv^2\right) - kp\right)$ $+ y(t) \dfrac{p \sqrt{2k\frac1{p} - 2k\frac1{y(t)} + v^2}}{2k + pv^2}$ (Wolfram|Alpha)
Now solving for $c_2$ in terms of $p$, $v$ and $k$ (at $t = 0$): $c_2 = \dfrac{k p^{3/2}}{(2k + pv^2)^{3/2}} \ln\left(p\sqrt{pv^2(2k + pv^2)} + p^2v^2 + kp\right) + \dfrac{p^2|v|}{2k + pv^2}$
Substituting back in, converting $\ln(a) - \ln(b)$ to $\ln(a/b)$: $\pm t = \dfrac{k p^{3/2}}{(2k + pv^2)^{3/2}} \ln\left(\dfrac{y(t) \left(\sqrt{(2k - 2k\frac{p}{y(t)} + pv^2)(2k + pv^2)} + 2k + pv^2\right) - kp}{p\sqrt{pv^2(2k + pv^2)} + p^2v^2 + kp}\right)$$+ \dfrac{y(t)p\sqrt{2k\frac1{p} - 2k\frac1{y(t)} + v^2} - p^2|v|}{2k + pv^2}$
And for $v < 0$, we need to switch the sign in front of $t$ on the left hand side (multiply by $\frac{|v|}{v}$) to preserve initial conditions: $t = \dfrac{k p^{3/2}}{(\frac{2k}{v^2} + p)^{3/2}v^3} \ln\left(\dfrac{y(t) \left(\sqrt{(2k - 2k\frac{p}{y(t)} + pv^2)(2k + pv^2)} + 2k + pv^2\right) - kp}{p\sqrt{pv^2(2k + pv^2)} + p^2v^2 + kp}\right)$$+ \dfrac{y(t)pv\sqrt{\frac{2k}{pv^2} - \frac{2k}{y(t)v^2} + 1} - p^2v}{2k + pv^2}$
So, after all this. The function becomes undefined when $2k < -pv^2$ because of the $\sqrt{2k + pv^2}$, but numeric integration clearly shows that there should be values just as legitimate here. What should I do to fix this? What's missing?
Resolution: Joriki's answer pointed me in the right direction. After wading through mounds of complex algebra, I have verified that the result above is correct for complex numbers and that to define the function in terms of real numbers when $2k < -pv^2$ you must use the arctangent function ($\ln(x + iy) = \ln(\sqrt{x^2 + y^2}) + i\tan^{-1}(\frac{y}{x})$, see Complex Logarithm). So here is my result for $2k < -pv^2$:
$t = -\dfrac{k p^{3/2}}{(-\frac{2k}{v^2} - p)^{3/2}v^3} \tan^{-1}\left(\dfrac{(k+pv^2)\sqrt{(-2k-pv^2)(-2k\frac{p}{y} + 2k + pv^2)} - (-k\frac{p}{y} + 2k + pv^2)\sqrt{pv^2(-2k - pv^2)})}{(k+pv^2)(-k\frac{p}{y} + 2k + pv^2) - (2k+pv^2)\sqrt{pv^2(-2k\frac{p}{y} + 2k + pv^2)}}\right)$$+ \dfrac{y(t)pv\sqrt{\frac{2k}{pv^2} - \frac{2k}{y(t)v^2} + 1} - p^2v}{2k + pv^2}$
(Any suggestions for simplification are appreciated!)
Edit: Fixed.