I have an ODE and the initial conditions for it. At those conditions, the function exists but the derivative approaches infinity. The equation for y(x), with K $\geq$ 0 and x$\geq$ 0,
$K=\frac{1}{y\sqrt{\dot{y}^2+1}}+\frac{\ddot{y}}
{(\dot{y}^2+1)^{3/2}}$,$y(0)=0, \dot{y}(0)=\infty$
This equation is the shape of a meniscus a pressure differential (Young-Laplace) for reference on its physical meaning. The initial conditions correspond to an axisymmetric shape around the x-axis. There is no analytical solution to this equation in most cases, so solving numerically is required. It can be split into two equations as:
$\dot{y}_1=y_2$,
$\dot{y}_2=K (y_2^2+1)^{3/2}-\frac{(y_2^2+1)}{y_1}$,
Attempting to solve this set of equations numerically using the initial conditions above is challenging, and nearly impossible. My machine precision $(\epsilon)$ is 2.22e-16, and using initial conditions of $y(0)=\epsilon, \dot{y}(0)=1/\epsilon$ provides a reasonable, yet incorrect answer. The solution is stopped when y(?)=1, if it ever gets that far before numerical errors take over. I've been using Matlab and tried the simple tricks like increasing tolerances, changing solvers (ode15s, 113, 23s, 45, 89) without luck.
In short my question is: Can this equation be rewritten in such a way that numerical solving is possible starting from the symmetry axis and/or it behaves nicely?
Hopefully something along these lines numerical solution of ode singularity that I'm not creative enough to see. Transforming from y(x) to x(y) makes the initial conditions reasonable, but the ODE is no better behaved and is still left with a singularity at y(0)=0.