# Stiffness and NDSolve error NDSolve::ndsz:

I have a 2nd Order ODE problem I am trying to solve numerically in Wolfram:
$$(k^2-1+\frac{2x}{x^2+c^2})\frac{d^2y}{dx^2}+\frac{c^2-x^2}{c^2+x^2}\frac{dy}{dx}+x\frac{x^2-3c^2}{(x^2+c^2)^3}y=0$$ where $$c=0.7$$,$$k=\sqrt{1-\frac{2b}{b^2+c^2}}$$ for $$b=50$$. My initial conditions are $$y(b)=1$$ $$\frac{dy}{dx}(b)=0$$.

I get the following error:

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

My code is here:

sola =
NDSolve[
{(k^2 - 1 + 2*x/(x^2 + c^2))*
D[y[x], x, x] + ((c^2 - x^2)/(x^2 + c^2)^2)*
D[y[x], x] + (x*(x^2 - 3 c^2)/(x^2 + c^2)^3)*y[x] == 0,
y[b] == 1, y'[b] == 0},
y, {x, 0, 100}]


Any idea of what to do here?

• At $x=b$, the coefficient of $y''$ vanishes; hence, you have a singularity. May 15, 2019 at 0:09

We can use the method of the false transient to solve this problem. We transform the equation to the form of the heat equation and then integrate it on both sides of the singular point using FEM:

Needs["NDSolveFEM"]

c = .7; b = 50; k = Sqrt[1 - 2*b/(b^2 + c^2)]; x0 = 100; x1 = 50;
sola = NDSolveValue[{-D[y[t, x], t] +
D[y[t, x], x,
x] + (((c^2 - x^2)/(x^2 + c^2)^2)*
D[y[t, x], x] + (x*(x^2 - 3 c^2)/(x^2 + c^2)^3)*
y[t, x])/(k^2 - 1 + 2*x/(x^2 + c^2)) ==
NeumannValue[0, x == x1], y[0, x] == 0,
DirichletCondition[y[t, x] == 1, x == x1]},
y, {t, 0, 10}, {x, x1, x0},
Method -> {"FiniteElement", "InterpolationOrder" -> {y -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.1}}];
solb = NDSolveValue[{-D[y[t, x], t] +
D[y[t, x], x,
x] + (((c^2 - x^2)/(x^2 + c^2)^2)*
D[y[t, x], x] + (x*(x^2 - 3 c^2)/(x^2 + c^2)^3)*
y[t, x])/(k^2 - 1 + 2*x/(x^2 + c^2)) ==
NeumannValue[0, x == x1], y[0, x] == 0,
DirichletCondition[y[t, x] == 1, x == x1]},
y, {t, 0, 10}, {x, 0, x1},
Method -> {"FiniteElement", "InterpolationOrder" -> {y -> 2},
"MeshOptions" -> {"MaxCellMeasure" -> 0.1}}];


The solution quickly converges on t and it looks like this

{Plot[Evaluate[Table[sola[t, x], {t, 1, 10, 1}]], {x, x1, x0},
PlotRange -> All],
Plot[Evaluate[Table[sola[t, x], {t, 1, 10, 1}]], {x, x1, x1 + 1},
PlotRange -> All]}

{Plot[Evaluate[Table[solb[t, x], {t, 1, 10, 1}]], {x, 0, x1},
PlotRange -> All],
Plot[Evaluate[Table[solb[t, x], {t, 1, 10, 1}]], {x, x1 - 1, x1},
PlotRange -> All]}


Combining solutions we find

Plot[If[x < 50, solb[10, x], sola[10, x]], {x, 0, 100},
PlotRange -> All]
`

• Could please recommend any material (Books, pdf...) where I may read about this method applied to Ordinary differential equations like this example? May 20, 2019 at 14:51
• @Ninpou Read this site first mathematica.stackexchange.com/… May 20, 2019 at 14:55