Can this gravitational field differential equation be solved, or does it not show what I intended?

Yes you're modeling it wrong. Gravity is an attractive force, so if the mass is at r = 0, and r is nonnegative, then the gravity force should point towards zero, i.e. negative.

$$ - \frac{GMm}{r^2} = m\ddot r$$

(This assumes $m\ll M$ so we don't need to consider reduced mass. See the other answers for the more general case.)

But now Alpha even refuses to solve it. That's fine, because Mathematica can't solve many things. The usual approach is to notice that

$$ \frac{d^2r}{dt^2} = \frac{dv}{dt} = \frac{dv}{dr}\frac{dr}{dt} = v \frac{dv}{dr} $$

to convert that 2nd-order ODE to a 1st-order ODE:

$$ -GM \frac{dr}{r^2} = v \,dv $$

This gives the solution v(r).

$$ \frac{GM}r + \text{constant} = \frac{v^2}2 \qquad (*) $$

But recall $v = dr/dt$. So we have another 1st-order ODE to solve, which gives r(t). The solution should look like that “hell of answer” because of the square root.

 

Note: If you rearrange the terms you should see that Equation (*) is just conservation of energy.

I think you are doing it all wrong. Your Newtonian gravitational force equation looks wrong to me in spherical coordinates.

Here's a different approach using Lagrangian and Hamiltonian approach. The solution is

$$\dfrac{\mathrm dr}{\mathrm dt}=\sqrt{\dfrac{2E}\mu-\dfrac{p_\phi^2}{\mu^2}\dfrac1{r^2}-\dfrac{2V(r)}\mu}$$

Where

$$ V(r) = -\frac1r $$

The correct differential equation is: $$ - \frac{G M m}{r^2} = \mu \ddot r $$ where $ \mu = \frac{M m}{M+m} $ is the reduced mass.

This can be simplified by dividing by $\mu$ $$ - \frac{G(M+m)}{r^2} =\ddot r $$

The problem can be further simplified by setting: $ G(M+m) $ equal to a constant, usually 1/2. $$ - \frac{1}{2} = r^2 \ddot r $$

This is a 2nd order quasilinear nonhomogenous ordinary differential equation. This equation is ill-formed, it does not have a unique solution. Adding initial or boundary conditions will lead to a unique particular solution.

Like the more general Kepler orbits, radial orbits can also be classified as elliptic, parabolic, or hyperbolic, corresponding to three forms of the particular solutions.

In the parabolic case, setting $ G(M+m) = 2/9 $, with initial conditions $ r(1)=1, \ \dot r(1)=2/3 \ $ leads to a simple solution:$$ t = r^{3/2} $$

In the elliptic case, setting $ G(M+m) = 1/2 $, with initial conditions $ r(\pi/2)=1, \ \dot r(\pi/2)=0 \ $, the particular solution is: $$ t = \arcsin(\sqrt{r})-\sqrt{r(1-r)} $$

In the hyperbolic case, setting $ G(M+m) = 1/2 $, with initial conditions $ t_0 = \sqrt{2}-\operatorname{arcsinh}(1)$, $ r(t_0)=1 $, $ \dot r(t_0)=\sqrt{2} $, the particular solution is:$$ t = \sqrt{r(1+r)} - \operatorname{arcsinh}(\sqrt{r}) $$