# 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? – J. M. will be back soon Mar 21 '18 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.) – Michael E2 Mar 21 '18 at 23:53
• @J.M. I think inside ParametricPlot, N[DSolve[..]] is tried, which now calls NDSolve, when DSolve fails. – Michael E2 Mar 21 '18 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. – bbgodfrey Mar 22 '18 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 Mar 22 '18 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 Mar 22 '18 at 22:13
• The problem seems to be the singularity x==0. If you set t1=1 everthing works! – Ulrich Neumann Mar 23 '18 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 . – Ulrich Neumann Mar 23 '18 at 7:22
• Thank you for your answer! – John Mar 23 '18 at 20:33