# Singularity problem of a 2nd order differential equation

I would greatly appreciate your help in implementation of Mathematica for solving a 2nd order differential equation expressing some important engineering problem.

Apparently, this equation has no analytical solution. Therefore I attempt to solve it numerically with maximal possible approach to the analytical effect, while using the ParametricNDSolve command.

The equation has 3 variables: n,Ω and ψ and differentiation is by z ($0<z<1$). n can be above zero and below it, however the values rarely exceed ±3. Ω and ψ are positive numbers and they can reach big values. As the 1st approach I have defined n=3. The expression is:

Solution=ParametricNDSolve[{Ω*f''[z][z]-f'[z]-ψ*f[z]^n==0,f-Ω*f'==1,f'==0},f,{z,0,1},{Ω,ψ}];


Unfortunately Mathematica is unable to solve this problem in wide ranges of Ω and ψ and produces the following messages:

ParametricNDSolve::ndsz: At z$349 == 0.9049990200792034, step size is effectively zero; singularity or stiff system suspected. >> ParametricNDSolve::ndsz: At z$343 == 0.8954517351902076, step size is effectively zero; singularity or stiff system suspected. >>
ParametricNDSolve::ndsz: At z\$343 == 0.8548633486473108, step size is effectively zero; singularity or stiff system suspected. >>
General::stop: Further output of ParametricNDSolve::ndsz will be suppressed during this calculation. >>


And the results at some points are non-physical.

Actually, the practically important solution is at z=1. I have tried that in the command: {z,1,1} and did not get these messages. However they have appeared while attempting to create the 3 dimensional plot:

Plot3D[Evaluate[f[Ω,ψ]/.Solution],{Ω,1/5,20},{ψ,0,7}]


And no progress has been achieved.

After reading answers to other questions on the command ParametricNDSolve I have tried to define a solution method:

Method→"StiffnessSwitching","ExtrapolationHandler"→{Indeterminate&,"WarningMessage"→False}


Unfortunately it did not result in any progress as well.

I would greatly appreciate your opinions if something else can be done in this case in order to overcome the singularity. Thank you a lot in advance.

• There are syntax errors in your code, such as \[Psi]*f[z]  instead of ψf[z], \[CapitalOmega]*f''[z] instead of Ωf"[z]  and f''[z] instead of f"[z]. The latter means that ''  this is twice Prim, but not the quote sign. After I corrected all this, your equation nicely solves. I propose to close this question as the off-topic. The problem is due to basic syntax errors. – Alexei Boulbitch Sep 28 '15 at 9:50
• Dear Alexei, thank you a lot for your comment. Unfortunately, the expression just has not been correctly reproduced by copying from my Mathematica code and I have edited my question. In fact there are no syntax errors there, and the procedure runs, yielding correct results at some points. Unfortunately, there is a problem of singularity at other points, where the results are non-physical. – Sitra Ahra Sep 28 '15 at 9:59
• @SitraAhra For copying & formatting code, you may find this this meta Q&A helpful. – Michael E2 Sep 28 '15 at 10:32
• Dear Michael, thank you a lot for this useful comment. I will implement your important suggestion in my next questions. – Sitra Ahra Sep 28 '15 at 10:44

Fast and dirty, try this:

 Clear[\[CapitalOmega], \[Psi], sol];
n = 3;
\[Psi] = 1;
lst = Flatten[
Table[sol =
NDSolve[{\[CapitalOmega]*f''[z] - f'[z] - \[Psi]*f[z]^n == 0,
f - \[CapitalOmega]*f' == 1, f' == 0}, f, {z, 0, 1}];
Table[
{\[CapitalOmega], z, f[z] /. sol}, {z, 0, 1,
0.1}], {\[CapitalOmega], 1, 2, 0.1}] /. {x_, y_, {z_}} -> {x, y,
z}, 1]


Now one may plot the obtained list, lst:

    ListPlot3D[lst,
AxesLabel -> {Style["\[CapitalOmega]", 18, Italic],
Style["z", 18, Italic], Style["Res", 18, Italic]}]


yielding this: Have fun!

• Dear Alexei, thank you a lot! Unfortunately this procedure does not solve the singularity problem. For example taking [Psi] = 4 (and bigger) the same error messages appear and some results are again non-physical. At the same time the problem could be simplified by considering z=1 only. The problem is that perhaps I do not know how to do it. Or perhaps they are some addtional methods. Thank you a lot again. – Sitra Ahra Sep 28 '15 at 13:39
• @Sitra Ahra Well, it seems that this equation, indeed, becomes singular at \[Psi]=3.92... Nothing to do. I tried the Method -> "StiffnessSwitching" ` it does not help. – Alexei Boulbitch Sep 28 '15 at 14:36