# Numerically solve 2nd order differential equation with singularity

Consider a second order differential equation with a potential that diverges at some generic value in the variable. For example:

$$-y^{\prime\prime}(s)+\frac1{\mathrm{cn}{(s\mid k^2)}}y(s)=0$$

where $\mathrm{cn}(s\mid k^2)$ is the JacobiCN[s, k^2] function. Now, suppose, the differential equation has to be solved on the interval of a full period, namely $s\in[0,2\mathbb{K}]$, where $\mathbb{K}$ is the complete elliptic integral of the first kind EllipticK[k^2], with the initial conditions $y(0)=0$ and $y^\prime(0)=1$.

If I just naïvely start the NDSolve function for, say, $k^2=0.7$:

NDSolve[{-y''[s] + y[s]/JacobiCN[s, 0.7] == 0, y[0] == 0, y'[0] == 1},
y[s], {s, 0, 2 EllipticK[0.7]}]


Mathematica rightfully complains that a singularity is detected at half of the evaluation range, namely at the value s = EllipticK[k^2] = 2.07536:

NDSolve::ndsz: At s == 2.0753631352632516, step size is effectively zero; singularity or stiff system suspected.

And also the resulting interpolating function terminates at this value:

{{y[s] -> InterpolatingFunction[{{0., 2.07536}},<>][s]}}
`

So, evidently, the solution does not overcome the point of singularity and is given only for the half of the required interval. I heard that there are numerical techniques to overcome singularities when solving differential equations. It makes me wonder that Mathematica does not activate these by default. What can I do to solve a differential equation beyond a singularity numerically in Mathematica?

-
This documentation about events might be helpful. –  Silvia May 9 '13 at 9:11
(I'm not posting an answer since this is rather specific to the DE concerned.) You have a coefficient that is a Jacobi elliptic function, which makes Picard's theory applicable (in complete analogy to the use of Floquet theory for differential equations with trigonometric coefficients, e.g. Mathieu). You might have better luck with an analytic approach instead of piggishly persisting with crossing through a singularity. –  Ｊ. Ｍ. May 9 '13 at 10:17
One other possibility you can do if you're unable or unwilling to follow through with Picard analysis is to do a bit of "pole-vaulting"; that is, integrate along a path in the complex plane that avoids the poles of $\mathrm{nc}(u\mid m)$. –  Ｊ. Ｍ. May 9 '13 at 12:36
"...involves a much more complicated, but still strictly elliptic potential term." - but you at least know where the poles might possibly be, yes? If so, your only job is to find a nice contour that avoids those nasty poles. In any event, if you want worked examples of Picard analysis, you might want to look up treatments of the well-studied Lamé equation. –  Ｊ. Ｍ. May 9 '13 at 13:15
I will just quote "A powerful way of computing the twelve Jacobian elliptic functions for real or complex values of both the argument z and the modulus k is to use the definitions in terms of theta functions given in ..." as found in NIST Handbook of Mathematical Functions pg. 566 –  Spawn1701D May 9 '13 at 13:34