Since people seem to be interested in the topic, I'll say a few more words.
1) As various people (including me) have said, the best approach is to
formulate the problem as a function of time: h(t), h'(t), and h''(t).
This is the conventional approach, and it has lots of advantages:
closer to the real physics, symplectic solver, et cetera.
Pedagogical point: It's always better to emphasize the good answer,
and explain why it is good, especially in the introductory course.
There will be plenty of time later for experts to delve into why the
bad answers are bad.
2) Having said that, let's put on our wizard hats and look into the
not-recommended version. It can be salvaged, and the salvage process
sheds considerable light on the root cause of the problem. Here's the
recipe:
The simplified equation of motion is:
h' = √(2 g h) [1]
The RHS is clearly a function of h. Arguably it is indirectly, implicitly
a function of t, since h itself depends on t, but that is useless to us
for present purposes, since we don't yet know the functional form of h(t).
The whole point of the exercise is to *find* h(t).
So there is nothing to lose and much to gain by thinking of it as a function
of h. We could even go so far as to write:
dt/dh = 1/√(2 g h) [2]
and then solve for t(h). In this case there is a convenient solution using
calculus (rather than numerical methods):
∫ dt = ∫ 1/√(2 g h) dh [3]
but let's pretend we didn't notice that.
Obviously integrating equation [2] or [3] numerically is a mess because the
RHS is infinite at h=0. However, we can salvage the situation as follows:
*Tabulate* dt/dh as a function of h, over the range of interest. The graph
is here: https://www.av8n.com/physics/img48/dt-dh-integrator.png
Our mission, should we decide to accept it, is to find the area under the curve.
(Or in this case, the shaded area to the left of the curve. I plotted h vertically.)
To find the area between h=1 and h=2 we could take the value of dt/dh at h=1
and integrate downward, or we could take the value at dt/dh at h=2 and integrate
upward, or (!) take the value at the midpoint of the interval (h=1.5), which is
a superior estimate of the integrand.
When finding the area between h=0 and h=1 we have only two of those three
options. The integrand is infinite at h=0. However the values at h=1 and
(especially!) h=0.5 are still viable estimates. It's an improper integral
because the integrand diverges, but it is integrable.
End of recipe.
Note that this trick of using an improved estimate of the integrand does not work
if you think of everything as a function of t. That would require knowing h(t)
or at least h'(t), which we don't know. Remember, finding h(t) is the whole point
of the exercise.
This recipe is not something I cooked up special to handle this case. It is
standard good practice to make a graph of what you know. This is mentioned
explicitly in Apostol's _Calculus_ book and probably also in Acton. In a
higher-dimensional space, it turns into a /flow field/ with little vectors
at each point showing which way the derivative goes.
Make a habit of doing this. Encourage students to do this. The harder the
problem, the more useful this is.
It's the same thing with electronics: Draw the circuit diagram. It never
ceases to amaze me when beginners try to figure things out in their head in
situations where world-class experts would draw the diagram.
===
Related point: If you think of it in terms of h(t) and h'(t) and time-stepping
starting from t=0, it is clearly a /stiff/ differential equation, because there
is no possible time-step during which h' does not change significantly during
that interval. It changes by a factor of infinity during that first step.
===
Additional point, affecting both the dh/dt and dt/dh versions: Whenever you
write down a square root, ask yourself why it's not the plus-or-minus square
root. In this case, the energy principle that was used to derive the equation
of motion does *not* distinguish between plus and minus velocity. And at the
critical point, both branches coincide at zero. You know that you are only
interested in positive velocities, but the computer can't read your mind.
Writing the positive square root gives away the answer at all h>0 but *not*
at h=0.
Rule of thumb: Beware the proximity of other solutions. Even if you know you
are not interested in them, the solver doesn't know that (unless you intervene
using very sophisticated methods). This is a fairly common way to stumble
into a stiff differential equation.
In the present case, starting at zero, the solver would be equally happy to
integrate backwards or forwards in time, at least for the first step. This
symmetry is one more way of explaining why it gets stuck.