# Specific problem with NDSolve step size, stiffness

I am trying to solve a non-linear second order boundary value problem in a finite interval. The differential equation is $$y''-\frac{a}{b}y-\frac{u_n}{b}y^3-\frac{ge_0}{b}x=0,$$ with $$x\in[0,L]$$, along with the boundary conditions, $$y'(0)-\frac{o}{b}y(0)=0,\qquad y'(L)+\frac{l}{b}y(L)=\frac{ge_0L}{b}.$$ In the code attached below, with the chosen parameters, the solution is regular and shows up until $$a=1.2$$. On using $$a=1.3$$ and beyond, I get the error "At x == 0.99505895226764, step size is effectively zero; singularity or stiff system suspected."

I am changing the input right after $$z$$ in Plot[Evaluate...

Working:

ClearAll;
Clear[Derivative]
fneeonnd[z_?NumericQ, a_?NumericQ, b_?NumericQ, o_?NumericQ,
l_?NumericQ, g_?NumericQ, L_?NumericQ, c_?NumericQ, u_?NumericQ,
un_?NumericQ, eo_?NumericQ] := (
ss = NDSolve[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x ==
0, y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)},
y, {x, 0, L}
, Method -> "StiffnessSwitching", WorkingPrecision -> 100];
y[z] /. ss)
Plot[Evaluate[
fneeonnd[z, 1.2, 0.05, 0.1, 0.1, 0.5, 1, 0.1, 0.1, 0.4, 1]], {z, 0,
1}, PlotStyle -> {Thickness[0.007]}, BaseStyle -> {FontSize -> 20},
Frame -> {True}, FrameLabel -> {"z", \[Eta]},
FrameStyle -> {{Black, Directive[Thick]}, {Black,
Directive[Thick]}, {Black, Directive[Thick]}, {Black,
Directive[Thick]}}]


NOT Working:

ClearAll;
Clear[Derivative]
fneeonnd[z_?NumericQ, a_?NumericQ, b_?NumericQ, o_?NumericQ,
l_?NumericQ, g_?NumericQ, L_?NumericQ, c_?NumericQ, u_?NumericQ,
un_?NumericQ, eo_?NumericQ] := (
ss = NDSolve[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x ==
0, y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)},
y, {x, 0, L}
, Method -> "StiffnessSwitching", WorkingPrecision -> 100];
y[z] /. ss)
Plot[Evaluate[
fneeonnd[z, 1.3, 0.05, 0.1, 0.1, 0.5, 1, 0.1, 0.1, 0.4, 1]], {z, 0,
1}, PlotStyle -> {Thickness[0.007]}, BaseStyle -> {FontSize -> 20},
Frame -> {True}, FrameLabel -> {"z", \[Eta]},
FrameStyle -> {{Black, Directive[Thick]}, {Black,
Directive[Thick]}, {Black, Directive[Thick]}, {Black,
Directive[Thick]}}]


• Look at your ODE. The main terms are: y'' = y^3. What does this mean? y'' has to to with the curvature, that is, if y>0, the larger y, the more it is upwards bent and the faster it grows. The same for y<0. Therefor you can solve only for values where y stays small. Commented Nov 6, 2020 at 21:23
• As far as I could see, y was actually getting smaller as I was increasing the input a. One can obtain the plots for a=0.5,0.75...1.2 to see this. Commented Nov 6, 2020 at 21:26
• As Daniel Huber says, curvature explodes at certain x. It shows to depend strongly on parameter b. Chose b == .006 and you can go up with a to a ==2 . Commented Nov 7, 2020 at 4:39

Often, Method -> "Shooting" with a good guess for the "StartingInitialConditions" is necessary to solve ODE boundary value problems. Here, we use a similar do-it-yourself shooting approach. As an example, use the parameters given for the "not working" case in the question:

{a, b, o, l, g, L, c, u, un, eo} =
{1.3, 0.05, 0.1, 0.1, 0.5, 1, 0.1, 0.1, 0.4, 1};

sp = ParametricNDSolveValue[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x == 0,
y'[0] - (o/b)*y[0] == 0, y[0] == y0}, {y[x], y'[L] + (l/b)*y[L] - (g*eo*L/b)},
{x, 0, L}, {y0}];

FindRoot[Last[sp[y00]], {y00, -.02}, Evaluated -> False]
First[sp[y00 /. %]];
Plot[%, {x, 0, L}, ImageSize -> Large, AxesLabel -> {x, y},
LabelStyle -> {15, Bold, Black}, PlotRange -> All]
(* {y00 -> -0.041088} *)


Of course, a reasonable guess is needed for y00 in FindRoot, here somewhere in the range {-.01, -.08}. A more ambitious case, a = 3, also works, yielding

(* {y00 -> -0.0163729} *)


It turns out that solving the linearized ODE yields a sufficiently good seeding for the standard Method -> "Shooting" to produce solutions of the nonlinear ODE for all values of a. The example below provides results for interger values of a between 0 and 10.

Clear[a];
Table[
tem = DSolveValue[{y''[x] - (a/b)*y[x] - (g*eo/b)*x == 0,
y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)}, y[0], {x, 0, L}];
NDSolveValue[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x == 0,
y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)}, y[x], {x, 0, L},
Method -> {"Shooting", "StartingInitialConditions" -> {y[0] == tem,
y'[0] == (o/b) y[0]}}], {a, 0, 10}];
Plot[Evaluate@%, {x, 0, L}, ImageSize -> Large, AxesLabel -> {x, y},
LabelStyle -> {15, Bold, Black}, PlotRange -> All,
PlotLegends -> Placed[Range[0, 10], {.3, .6}]]


One can use the FEM method or use the FEM method to find good starting initial conditions. I do the latter in fneeonnd but save the FEM solution in solFEM. The check at the end shows that the shooting method did a better job of satisfying the boundary conditions.

fneeonnd[z_, a_?NumericQ, b_?NumericQ, o_?NumericQ, l_?NumericQ,
g_?NumericQ, L_?NumericQ, c_?NumericQ, u_?NumericQ, un_?NumericQ,
eo_?NumericQ] := (
Clear[x, y];
solFEM =
NDSolveValue[{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x ==
NeumannValue[(o/b)*y[x], x == 0] +
NeumannValue[-(g*eo*L/b) + (l/b)*y[x], x == L]},
y, {x} \[Element] Line[{{0}, {L}}]];
NDSolveValue[
{y''[x] - (a/b)*y[x] - (un/b)*(y[x])^3 - (g*eo/b)*x == 0,
bcs = {y'[0] - (o/b)*y[0] == 0, y'[L] + (l/b)*y[L] == (g*eo*L/b)}},
y, {x, 0, L},
Method -> {"Shooting",
"StartingInitialConditions" -> {y[0] == solFEM[0],
y'[0] == solFEM'[0]}}(*,WorkingPrecision\[Rule]100*)])

solFN = fneeonnd[z, 1.3, 0.05, 0.1, 0.1, 0.5, 1, 0.1, 0.1, 0.4, 1];
Plot[solFN[z], {z, 0, 1}, PlotStyle -> {Thickness[0.007]},
BaseStyle -> {FontSize -> 20}, Frame -> {True},
FrameLabel -> {"z", \[Eta]},
FrameStyle -> {{Black, Directive[Thick]}, {Black,
Directive[Thick]}, {Black, Directive[Thick]}, {Black,
Directive[Thick]}}]


Plot[{solFN[z], solFEM[z]}, {z, 0, 1}]


bcs /. Equal -> Subtract /. {{y -> solFN}, {y -> solFEM}}
(*  {{0., 2.55954*10^-6},        <-- Shooting
{-0.00153343, -0.0726405}}  <-- FEM      *)
`