# Error when using NDSolve for $\epsilon y'' - y' + y^2 = 1$

Error when using NDSolve for $$\epsilon y'' - y' + y^2 = 1$$ with $$0 and $$y(0) = \frac{1}{3}$$, $$y(1)=1$$

My attempt:

Subscript[ϵ, 1] = 0.1;
Subscript[eqn, 1] =
NDSolve[{Subscript[ϵ, 1]*y''[x] - y'[x] + y[x]^2 == 1,
y[0] == 1/3, y[1] == 1}, y, {x, 0, 1}];


Error:

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

However I see from this post https://math.stackexchange.com/a/3217456/697826 that it was numerically solved. What is wrong with my syntax?

Thank you to both of the submitted Answers. They are equally helpful!

Nothing wrong with your syntax. This is because

1. The default setting of "Shooting" method doesn't work well on this problem.

2. There's an arguable backslide of "Shooting" method around v11. See this and this post for more info.

If you're in v12.0 or higher, in which the nonlinear FEM solver is added, the simplest work-around is to turn to FiniteElement:

Subscript[ϵ, 1] = 0.1;
solfem = NDSolveValue[{Subscript[ϵ, 1]*y''[x] - y'[x] + y[x]^2 == 1,
y[0] == 1/3, y[1] == 1}, y, {x, 0, 1}, Method -> FiniteElement]

Plot[solfem[x], {x, 0, 1}]


If you're before v11, then the following should work:

Subscript[ϵ, 1] = 0.1;
solshoot = NDSolveValue[{Subscript[ϵ, 1]*y''[x] - y'[x] + y[x]^2 == 1,
y[0] == 1/3, y[1] == 1}, y, {x, 0, 1},
Method -> {"Shooting",
"StartingInitialConditions" -> {y'[1] == 1, y[1] == 1}}]

solshoot // ListPlot


Notice I've made use of an undocumented syntax of ListPlot here. See this post for more info.

First, there are multiple solutions. Second, the one shown on math.SE has such a narrow range of feasible starting initial conditions fo $$\epsilon = 1/10$$ that the default shooting method in NDSolve is unlikely to find it, however the options are tweaked. It's probably best to manually implement the shooting method.

\[Epsilon] = 1/10;
psol = ParametricNDSolveValue[{\[Epsilon]*y''[x] - y'[x] + y[x]^2 ==
1, y[0] == 1/3, y'[0] == yp}, y, {x, 0, 1}, {yp}]

ypsol = FindRoot[psol[yp][1] == 1, {yp, -0.94, -0.93},
Method -> "Brent"]

(*  {yp -> -0.93371}  *)

sol = psol[yp] /. ypsol
ListLinePlot@sol


Here's another form of the shooting method, which allows one to figure out the interval to pass to the Brent method above:

Manipulate[
ListLinePlot@
NDSolveValue[{\[Epsilon]*y''[x] - y'[x] + y[x]^2 == 1, y[0] == 1/3,
y'[0] == yp}, y, {x, 0, 1}],
{yp, -0.95, -0.5, Appearance -> "Labeled"}]

• Thank you Michael! Commented May 3, 2021 at 4:17
• You're welcome, @Nalt! :) Commented May 3, 2021 at 12:06

Things are easy to solve, if you differentiate one time.

e1 = 1/10;
eqs0 = {e1*y''[x] - y'[x] + y[x]^2 == 1, y[0] == 1/3, y[1] == 1}
deq1 = D[eqs0[[1]], x]

ic2 = First@First@Solve[eqs0[[1]] /. x -> 1, y'[1]] /. y[1] -> 1 /.
Rule -> Equal

ysol = y /.
First@NDSolve[{deq1, ic2, y[0] == 1/3, y[1] == 1}, y, {x, 0, 1}]

{ysol[0], ysol[1]}
(*   {0.333333, 1.00001}   *)

Plot[ysol[x], {x, 0, 1}]
`