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?
