A differential equation Think of $t$ and $r$ as two independent variables. 


*

*Suppose $E$ be a function of $r$ and $V~$ be a function of $(t,r)$ such that both go to $0$ at $r=0$. 

*There exists a positive function $M(r)$ such that $M(0)=0$ and $V(t,r) = -\dfrac{M(r)}{R(t,r)}$ where $R$ is another positive function such that $R(0,r)=r$. 

*Let $p(r) = \dfrac{E(r)}{V(0,r)}$ be a function regular at $r=0$ such that $p(0) \in (-\infty,1)$. 

*Also define a function $a$ of $r$ such that, $a(r) = \dfrac{M(r)}{\dfrac{4}{3}\pi r^3}$. 
Then $a$ is also a positive definite function with a well-defined value at $r=0$.  

*Define $\alpha = a(0)$ 
Now look at this differential equation,
$$\frac{\dot{R}^2}{2} + V(t,r) = E(r)$$
Apparently this differential equation has a solution of the form,
$$\frac{t}{t_0} = \sqrt{\frac{\alpha}{a(r)}}\frac{F(p(r))}{F(p(0))} \left [1 - \left ( \dfrac{R(t,r)}{r} \right)^{\dfrac{3}{2}}~\cdot~\dfrac{F\left(~~ \dfrac{p(r)R(t,r)}{r} \right) }{F(p(r))} \right ] $$
where $t_0 = \sqrt {\dfrac{3}{8\pi \alpha}} F(p(0))$ 
and the function $F$ is defined over the interval $(-\infty,1)$ as,
$$F(x) = \left\{ 

\begin{array}{c c}   
-\frac{\sqrt{1-x}}{x} - \frac{1}{(-x)^{\frac{3}{2}}} \tanh^{-1} \left [ \sqrt{\frac{x}{x-1}} \right ] & x&lt;0 \\     
\frac{2}{3}  &  x =0 \\
\frac{1}{x^{\frac{3}{2}}}tan^{-1} \left [ \sqrt{\frac{x}{1-x}} \right ] - \frac{\sqrt{1-x}}{x}  & 0&lt;x&lt;1 
\end {array} 
\right. $$
How does one get the above solution? 
 A: I'm sort of just guessing here (partly based on the solution already found). By explicitly plugging in $V(t,r) = - \frac{M(r)}{R(t,r)}$, you arrive at the ordinary differential equation (for each fixed $r$) for $R$ as
$$ \dot{R}^2 - \frac{2M}{R} = 2 E $$
where $M$ and $E$ are constants in time. Now, re-scale the original equation by $t = \lambda s$ and $R = \mu \rho$. Then $\partial_t R = \frac{\mu}{\lambda} \partial_s \rho$. Then you can solve $(\frac{\mu}{\lambda})^2 = \frac{2M}{\mu} = 2E$ to reduce the equation to 
$$ \dot{\rho}^2 - \frac{1}{\rho} = 1 $$
(the weights $\mu$ and $\lambda$ will, roughly speaking, give you the weights $p(r)$ and $a(r)$ in your question). Now note that this scaling degenerates if $E = 0$. In the case that $E = 0$, the equation can be solved by quadrature:
$$ \dot{x}^2 = x^{-1} \Rightarrow \sqrt{x} dx = dt \Rightarrow x^{3/2} \sim t $$
This gives the solution to the homogeneous case. In the inhomogeneous case, you take that as a sort of integrating factor: assume that $\rho^{3/2} f(\rho) \sim t$, this implies that 
$$ \dot{\rho} \left( \rho^{3/2} f(\rho) \right)' = 1 $$
we plug this into the equation, which we first re-arrange as
$$ \dot{\rho} = \sqrt{\frac{1}{\frac{\rho}{1+\rho}}} $$
and we conclude that 
$$ \sqrt{\frac{\rho}{1+\rho}} = \frac{d}{d\rho}( \rho^{3/2} f(\rho)) $$
and you solve this by directly integrating it. I think this $f$ you find should be exactly the $F$ you wrote down above (I didn't check it myself). Note that due to the singular weight $\rho$ which degenerates as $\rho \to 0$, you will have to separately integrate in the regime where $\rho > 0$ and $\rho < 0$. The existence and uniqueness of solution is guaranteed by the theory of Fuchsian ODEs. 
