# 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. – Michael E2 May 15 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? – Ninpou May 20 at 14:51
• @Ninpou Read this site first mathematica.stackexchange.com/… – Alex Trounev May 20 at 14:55