I have the Cauchy problem:
$$
\frac{dx}{dt} = \frac{f(t, x)}{g(x) - t^3} \;,\qquad
x(t_0) = x_0 \;.
$$
It can be shown that when the solution reaches a vicinity of a certain point $(t_\ast, x_\ast)$ the r.h.s. behaves such a way:
$$
\frac{dx}{dt} \sim - \frac{\operatorname{sign}(t - t_\ast)}{|t - t_\ast|^{1/2}} \;.
$$
This singularity is integrable, so the solution looks like
I expect a series of such suingularities in the solution.
The problem is I can't integrate the ODE numerically. Neither explicit (Adams, RK, etc.) nor implicit (BDF) method can pass the singularity. Providing the jacobian doesn't help.
I tried to reformulate the problem, e.g. to define $x(t)$ as a parametric curve on the $(t, x)$ plane, $$ t = t(s) \;,\qquad y(s) = x(t) \;,\qquad (dt/ds)^2 + (dy/ds)^2 \equiv 1 \;. $$ The idea was that $s$ is the natural parametrization so the velocity along the curve is always constant. Doesn't help. The solver can not cope with the fact that $|dt/ds|$ and $|dy/ds|$ may differ from each other by many orders of magnitude.
Reading NumRecipes and such I can't find the proper method also.
UPD: Last news after some break. I couldn't find or construct any numerical method to resolve the caustic. The trick of @LutzL didn't help also (btw, the numerator and the denominator don't turn to zero simultaneously at any time, including caustic).
This problem arised as a self-similar model of a certain gas-dynamic problem. We've made a numerical gas-dynamical solver in a time-space domain and it turned out that a shock wave appears in the soution exactly at the same place the caustic appears (here $t$ is a self-similar variable). This is a kind of a happy end for us. However, the pure algoritmic issue still exists in the self-similar formulation of the problem. I don't know how to solve an ODE with a finite singularity at the r.h.s.
Thanks to everyone who cared.