# BEYOND Singularity with NDSolve

Here is the code I want to numerically integrate my nonlinear ODE"

sol = NDSolve[{y''[x] + y[x] == -1/(1 + y[x])^2, y[0] == 1, y'[0] == 1}, y, {x, 0, 5}]
Plot[y[x] /. sol, {x, 0, 3}]


NDSolve surely fails since it passes through the singular at the point x=2.532686648521246.

I tried another method from here to bypassing a singularity:

odey = R''[y] + R[y] + 1/(1 + R[y])^2;
BCm = {R[0] == 1, R'[0] == 1};
projODE = {odey /. R -> (p[#]/q[#] &) // Together // Numerator,
D[p[y]^2 + q[y]^2, {y, 2}]} == 0 // Thread;

projICS = BCm /. {R[y0_] == r0_ :> {p[y0], q[y0]} == {Numerator[r0], Denominator[r0]},
R'[y0_] ==
rp0_ :> ({D[p[y]/q[y], y] == rp0, D[p[y]^2 + q[y]^2, y] == 0} /.
y -> y0)};

rules = {AccuracyGoal -> Infinity, PrecisionGoal -> 20,
WorkingPrecision -> 50, MaxSteps -> 20000};
projSOL = NDSolve[{projODE, projICS}, {p, q}, {y, 0, 5}, rules,
Method -> {"EquationSimplification" -> {Automatic,
"SimplifySystem" -> True}}]

Plot[{p[y]/q[y] /. First@projSOL}, {y, 0, 4}, PlotLegends -> {p/q}]


but dosen't work.

My question is that, is there a convenient numerical continuation method to bypass the singularity so that NDSolve does not fail from range: {x,0,5} ?

• I think a good place to start is to try to copy/paste the code Michael posted in the thread you refer to. I tried and ran into a mess: missing variables and lots of errors. So I think first get his code working using the problem in that thread, then apply it to yours. For example, variable "rat" is not defined anywhere in that thread but I think it can be 10^-10 for example.
– josh
Dec 16, 2021 at 23:11

You need a method that is adapted to stiff DE. E.g. "LinearlyImplicitEuler"

sol = NDSolve[{y''[x] + y[x] == -1/(1 + y[x])^2, y[0] == 1,
y'[0] == 1}, y, {x, 0, 5}, Method -> "LinearlyImplicitEuler"]
Plot[y[x] /. sol, {x, 0, 3}, PlotRange -> All]


• Clever solution, fits well to the phaseplot in my answer! Dec 17, 2021 at 11:24
• The result of Plot[y[x] /. sol, {x, 0, 5}, PlotRange -> All] looks very strange. Dec 17, 2021 at 11:56
• Could you also look at what happens if you try Method -> "StiffnessSwitching"? Dec 17, 2021 at 15:19

To long for a comment:

Perhaps StreamPlot gives an idea of the possible solution?

StreamPlot[{yp, -y - 1/(1 + y)^2}, {y, -2, 1}, {yp, -3, 3},FrameLabel -> {y[x], y'[x]}]