# At t == ..., step size is effectively zero; \ singularity or stiff system suspected

I need to get a phase portrait for a non-linear oscillator, for this I wrote down the corresponding equations. 3 equations for different "a" which can take values from minus to plus infinity. In each equation, the last term is an expression for the nonlinearity of our oscillator. In general, the problem (error message NDSolve::ndsz) arises when I try to solve equations for "a" which value is more than zero. For example, if i change last term for another like that: (a*x)/(a+x), it works. Can't anderstand where is my mistake in the first case? My code is:

δ = 0.1;
ω0 = 1;
α1 = 0.1;
α2 = 0;
α3 = -0.5;
Subscript[t, 0] = 0;
Subscript[t, 1] = 100;

x1 = DSolve[{x''[t] +
2δ*x'[t] + (ω0^2)*x[t] + α1/x[t]^2 == 0,
x[0] == 1, x'[0] == 0}, x, {t, Subscript[t, 0], Subscript[t, 1]}];

x2 = DSolve[{x''[t] +
2δ*x'[t] + (ω0^2)*x[t] + α2/x[t]^2 == 0,
x[0] == 1, x'[0] == 0}, x, {t, Subscript[t, 0], Subscript[t, 1]}];

x3 = DSolve[{x''[t] +
2δ*x'[t] + (ω0^2)*x[t] + α3/x[t]^2 == 0,
x[0] == 1, x'[0] == 0}, x, {t, Subscript[t, 0], Subscript[t, 1]}];

Show[ParametricPlot[{x[t], x'[t]} /. x1, {t, 0, 60},
PlotRange -> {{-1, 1}, {-2, 2}}, PlotStyle -> Blue,
PlotLegends -> {"α = 0.1"}],
ParametricPlot[{x[t], x'[t]} /. x2, {t, 0, 60},
PlotRange -> {{-1, 1}, {-2, 2}}, PlotStyle -> Black,
PlotLegends -> {"α = 0"}],
ParametricPlot[{x[t], x'[t]} /. x3, {t, 0, 60},
PlotRange -> {{-1, 1}, {-2, 2}}, PlotStyle -> Red,
PlotLegends -> {"α = -0.5"}]]

• Only NDSolve[] and ilk can produce the NDSolve::ndsz error, yet you are using DSolve[] in your code. What am I missing? Commented Mar 21, 2018 at 23:52
• You might try tracking down which command give the error and restrict your code to what is necessary. If DSolve doesn't work, try NDSolve, if a numeric solution would be acceptable. (It seems x1 and x3 fail.) Commented Mar 21, 2018 at 23:53
• @J.M. I think inside ParametricPlot, N[DSolve[..]] is tried, which now calls NDSolve, when DSolve fails. Commented Mar 21, 2018 at 23:54
• Yes, DSolve is unable to solve the x1 and x3 ODEs. Incidentally, because the ODEs are autonomous, they can be reduced to first order, although it is unclear whether that would help. Also, I would guess that, if NDSolve eventually is called, it would fail at the turning points of the nonlinear oscillators or when x[t] == 0. Commented Mar 22, 2018 at 2:15
• I'm sorry, I used NDsolve and it doesn't work too, than I tried to use Dsolve to see what changes and forgot to change it back
– John
Commented Mar 22, 2018 at 8:31

Just use ParametricNDSolveValue to solve your problem:

sol = ParametricNDSolveValue[{x''[t] +2 \[Delta]*x'[t] + \[Omega]0^2*x[t] +\[Alpha] x[t]^2 == 0,x[0] == 1, x'[0] == 0}, x, {t, t0, t1}, {\[Alpha]}]
(*ParametricFunction[ <> ] *)

ParametricPlot[Evaluate[Table[{sol[\[Alpha]][t],
sol[\[Alpha]]'[t]}, {\[Alpha], {\[Alpha]1, \[Alpha]2, \[Alpha]3}}]], {t, t0,t1}, PlotRange -> All]


• Everything works, but you wrote last term as α*x[t]^2 but i need α/x[t]^2 and if change last it doesn't work.
– John
Commented Mar 22, 2018 at 22:13
• The problem seems to be the singularity x==0. If you set t1=1 everthing works! Commented Mar 23, 2018 at 7:06
• Check the equilibrium of your ode 0 == \[Omega]0^2*x0 + \[Alpha] /x0^2: Only for \[Alpha] <0 a real solution x0>0exists . Commented Mar 23, 2018 at 7:22