Analytic solution to Newtonian gravity differential equation

There actually exists 2 solutions.

We first make a change of variable $t=\log(s)$ to eliminate the troublesome -Infinity, here I'll use DChange for the task:

eq = {z''[t] + 1/z[t]^2 == 0, z'[-inf] == 0, z[0] == 0};

neweq = DChange[eq, t == Log@s, t, s, z[t]] /. inf -> Infinity
(* {1/z[s]^2 + s (z'[s] + s z''[s]) == 0, 
    z'[0] == 0, z[1] == 0} *)

general = DSolve[neweq[[1]], z[s], s]
(*{Solve[……], Solve[……]}*)

DSolve gives 2 unsolved equations as the general solution. By substituting the b.c. $z(1)=0$ into the equation:

soleq = general[[All, 1]];
soleq /. s -> 1 /. z[1] -> 0
(* {0 == C[2], 0 == C[2]} *)

C[2] is determined. The next task is to make use of the b.c. $z'(0)=0$. We differentiate the equation and solve for $z'(s)$:

First@Solve[D[#, s], z'[s]] & /@ (soleq /. C[2] -> 0)
(* {{z'[s] -> (Sqrt[2] Sqrt[1 + C[1] z[s]])/(s Sqrt[z[s]])}, 
   {z'[s] -> -((Sqrt[2] Sqrt[1 + C[1] z[s]])/(s Sqrt[z[s]]))}} *)

Given $z'(0)=0$, 1 + C[1] z[s] /. s -> 0 can only be equal to 0 i.e. z[s] -> -1/C[1] /. s -> 0. Substitute this condition back into soleq and solve for C[1]:

soleq /. C[2] -> 0 /. z[s] -> -1/C[1]
(* {(I π)/(2 Sqrt[2] C[1]^(3/2)) + Log[s] == 0, 
    -((I π)/(2 Sqrt[2] C[1]^(3/2))) + Log[s] == 0} *)
First@Solve[#, C[1]] & /@ % /. s -> 0
(* {{C[1] -> 0}, {C[1] -> 0}} *)

Substitute the value of C[1] and C[2] into the equation (for C[1] we need to take the limit), solve for z[s] and change the variable back:

0 == Limit[Subtract @@ #, C[1] -> 0] & /@ soleq /. C[2] -> 0
(* {0 == Log[s] - 1/3 Sqrt[2] z[s]^(3/2), 
    0 == Log[s] + 1/3 Sqrt[2] z[s]^(3/2)} *)
newsol = First@Solve[#, z[s]] & /@ %
(* {{z[s] -> (3^(2/3) Log[s]^(2/3))/2^(1/3)}, 
    {z[s] -> (3^(2/3) (-Log[s])^(2/3))/2^(1/3)}} *)

{sol1[t_], sol2[t_]} = z[s] /. newsol /. Log@s -> t
(* {(3^(2/3) t^(2/3))/2^(1/3), (3^(2/3) (-t)^(2/3))/2^(1/3)} *)

Both of the solutions satisfy the equation and the b.c.s:

eq /. {{z -> sol1}, {z -> sol2}} /. inf -> Infinity
(* {{True, True, True}, {True, True, True}} *)

It seems that the differential equation solver gets hung up on the boundary conditions for some reason. Let us therefore do some intermediate steps by hand. First, consider the general solution of the differential equation:

DSolve[{z''[t] + 1/z[t]^2 == 0}, z[t], t]

Solve[(-(Log[1 + (C[1] + Sqrt[C[1]] Sqrt[C[1] + 2/z[t]]) z[t]]/C[1]^(3/2)) + (Sqrt[C[1] + 2/z[t]] z[t])/C[1])^2 == (t + C[2])^2, z[t]]

This result means that an ordinary equation must be satisfied by z[t] in order to give the differential equation solution. Let's rewrite it as expr==0 with

expr = (-(Log[1 + (C[1] + Sqrt[C[1]] Sqrt[C[1] + 2/z[t]]) z[t]]/C[1]^(3/2))
        + (Sqrt[C[1] + 2/z[t]] z[t])/C[1])^2 - (t + C[2])^2;

Note that t only enters as t^2, which means that we can immediately replace t by its absolute value Abs[t] in all further considerations, since time in physics is a real variable (For simplicity, let us concentrate on the region where Abs[t]=t for now, and analytically continue the solution to the other region when we are done). We also see that there are two integration constants C[1],C[2] that we should use to obtain the boundary conditions we are interested in. First condition for t=0:

Series[expr /. t -> 0, {z[0], 0, 0}] // Normal

-C[2]^2

implies C[2] -> 0. For the second condition we need the first derivative z'[t]:

zp = z'[t] /. Solve[D[expr /. C[2] -> 0, t] == 0, z'[t]][[1]] // Simplify

(t C[1]^(3/2) Sqrt[ C[1] + 2/z[t]])/(-Log[ 1 + (C[1] + Sqrt[C[1]] Sqrt[C[1] + 2/z[t]]) z[t]] + Sqrt[C[1]] Sqrt[C[1] + 2/z[t]] z[t])

This almost looks like C[1] -> 0 should be chosen. However, we should take this limit carefully:

zpBC = Series[zp, {C[1], 0, 0}] // Normal

(3 t)/z[t]^2

We see that the first derivative is non-zero for C[1] -> 0, but we still need to take t -> -Infinity. It makes sense that z[t] should go to infinity in that limit as well (particle falling in from far away), but we need to make sure that it approaches infinity quicker than Sqrt[t] for the above first derivative to vanish. Indeed, we find:

Series[expr /. C[2] -> 0, {C[1], 0, 0}] // Normal

-t^2 + (2 z[t]^3)/9

So that

z[t] == (3^(2/3) t^(2/3))/2^(1/3)

And the first derivative properly goes to zero as t -> -Infinity:

zpBC /. z[t] -> (3^(2/3) t^(2/3))/2^(1/3)

2^(2/3)/(3^(1/3) t^(1/3))

Therefore, z[t] == (3^(2/3) t^(2/3))/2^(1/3) is likely the solution you are after. We can also check explicitly:

z[t_] := (3^(2/3) t^(2/3))/2^(1/3)
z''[t] + 1/z[t]^2

0

that the differential equation is indeed satisfied. Also, recall that by t we actually mean Abs[t]. This implies that substituting t -> -t also gives a solution. This one is actually the one you need, since for t<0 we have Abs[t] == -t and you expect the solution to be real on physical grounds.

Here is a very short answer, using energy conservation:

energy = 1/2 z'[t]^2 - 1/z[t];

zSolution[t_] = 
 z[t] /. Simplify@First[DSolve[energy == 0 && z[0] == 0, z[t], t]]

(* ==> (3^(2/3) (-t)^(2/3))/2^(1/3) *)

Energy conservation reduces the order of the differential equation. The energy at $t\to-\infty$ must be equal to zero if the velocity is zero in that limit. This follows from the fact that the limit $t\to-\infty$ can only be defined provided that the motion is non-periodic in time. This implies that the potential energy at $t\to-\infty$ must also vanish.

Added note about the energy:

To justify that the quantity called energy above is in fact conserved and therefore must obey the differential equation I use here, you could use the VariationalMethods package, even though it's overkill for this simple problem. The logic is: first I define lagrangian and show that its Euler-Lagrange equation of motion is in fact the desired differential equation. Having therefore described the problem equivalently by a Lagrangian, I can derive conservation laws from it. In particular, the conservation law I need is that of the energy, which in Mathematica is obtained as FirstIntegral[t] from the command FirstIntegrals:

Needs["VariationalMethods`"]

lagrangian = 1/2 Derivative[1][z][t]^2 + 1/z[t];

Simplify[EulerEquations[lagrangian, z[t], t]]    
(* ==> 1/z[t]^2 + (z'')[t] == 0 *)

energy = 
 FirstIntegral[t] /. FirstIntegrals[lagrangian, z[t], t]
(* ==> -(1/z[t]) + 1/2 z'[t]^2 *)

As you can see, this is the quantity I start with in the original post. So this is just a Mathematica-based way of deriving energy conservation, provided you're familiar with Lagrangian mechanics. Of course in real physics problems, you should usually identify the relevant conservation law before doing any computations at all. It saves a lot of work.