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

This is the equation I'm having trouble with:

$$G \frac{M m}{r^2} = m \frac{d^2 r}{dt^2}$$

That's the non-vector form of the universal law of gravitation on the left and Newton's second law of motion on the right. I assume that upon correctly modeling and solving this, I will have a function of time that gives the distance from a spherical mass in space (e.g. distance from the Earth from an initial condition of $r(0) = 10,000 \mathrm{km}$).

However, WolframAlpha gives a hell of an answer, which leads me to believe that I'm modeling this equation completely wrong. Can anyone shed some light on this problem?

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.

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}) $$

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 $$