# Solving ODE system with NDSolve: “singularity or stiff system suspected” [closed]

I'm trying to solve a set of 4 ODEs as specified below, however I'm getting this error

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


Free parameters:

m = 0; q = 1.3; α = 0.05; a0 = ((α^2 + 5 α + 6)/2)^0.5/q; c0 = 0.56;
ϵ = 0.0001;


Note ϵ is a small number used in place of 0, so can be changed and in fact the location of this singularity/ stiffness is dependent on ϵ. These are the four equations:

Clear[eqns, eq1, eq2, eq3, eq4 , sol]
Clear[a, b, c, g]

eq1 = a''[r] + (g'[r]/g[r] - c'[r]/2 + 2/r) a'[r] + (q^2 b[r]^2 Exp[c[r]])/g[r]^2 a[r] - m^2/g[r] a[r] == 0 ;
eq2 = b''[r] + (c'[r]/2 + 2/r) b'[r] - (2 q^2 a[r]^2)/g[r]^2 b[r] == 0;
eq3 = c'[r] + r a'[r]^2 + (r q^2 b[r]^2 a[r]^2 Exp[c[r]])/g[r]^2 == 0;
eq4 = g'[r] + (1/r - c'[r]/2) g[r] + (r b'[r]^2 Exp[c[r]])/4 - 3 r + (r m^2 a[r]^2)/2 == 0;


These are the boundary/ initial conditions:

eqns = {eq1, eq2, eq3, eq4,
a[ϵ] == a0, a'[ϵ] == ϵ,
b[ϵ] == ϵ , b'[ϵ] == ϵ,
c[ϵ] == c0, g[ϵ] == ϵ
};

sol = NDSolve[eqns, {a, b, c, g}, {r, ϵ, 2}]


The solution should qualitatively look like this

I'm not able to figure out where or how this stiffness shows up. Any sort of help would be appreciated, thank you!

EDIT: Originally used {φ, ϕ, χ, g} are now used as {a, b, c, g}

## closed as off-topic by István Zachar, MarcoB, user9660, xzczd, Bob HanlonDec 23 '15 at 17:56

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question cannot be answered without additional information. Questions on problems in code must describe the specific problem and include valid code to reproduce it. Any data used for programming examples should be embedded in the question or code to generate the (fake) data must be included." – István Zachar, MarcoB, Community, xzczd, Bob Hanlon
If this question can be reworded to fit the rules in the help center, please edit the question.

• I'd posted this question yesterday but after posting I found numerous mistakes in that so preferred to re-post a fresh question. – Skylar15 Dec 22 '15 at 5:30
• Under V10.1 with very slight changes to your code that shouldn't make any difference I do not get your error, g(r) looks parabolic or exponential going up to about 2 and with no "knee" while phi(r) "goes almost straight up" and is then horizontal at about 1.3601 starting about r=0.2 – Bill Dec 22 '15 at 6:39
• I've 10.0.2 running on OS X, can you please write what changes you did? Thanks! – Skylar15 Dec 22 '15 at 6:43
• I obtain the same sort of solution as @Bill with no change to the code. I use 10.0.3 on Windows 10. By the way, he probably did not see your comment to him, because you did not address it to him. – bbgodfrey Dec 22 '15 at 11:36
• @Bill and @bbgodfrey it's strange i'm getting that error and in fact NDSolve is giving me solution in the range {0.0001, 0.00022} only! – Skylar15 Dec 22 '15 at 11:41

First, this layer after layer of mandatory manual desktop-publishing followed by manual un-desktop-publishing and screwing around just to get greek characters on the screen and make code acceptably pretty adds layers of uncertainty and potential errors. For example, I suspect that somewhere in the editing process "curly phi" and "phi" got reversed in your post. Next, I would not be at all surprised that your actual real un-desktop-published code just might have a single incorrect greek character manual translation and that is why you are seeing one thing and everyone else is seeing something completely different. But we can't go back and kill Steve Jobs and there is nothing that will change this now.

Fresh start of MMA 10.1 and then I scrape-n-paste your code from your post with very minimal changes, other than all the required manual desktop publishing edits. This gives zero error/warning messages and the plots are what they are. Please scrape-n-paste-n-execute this code with zero manual edits to verify on your system and then track down the likely almost invisible difference between this and your code.

Clear[φ, ϕ, χ, g];
m = 0; qφ = 1.3; q = qφ/2; α = 0.05; χ0 = 0.56;
ϕ0 = Sqrt[(α^2 + 5 α + 6)/2]/qφ; ϵ = 0.0001;
eq1 = φ''[r] + (g'[r]/g[r] - χ'[r]/2 + 2/r) φ'[r] +
(qφ^2 ϕ[r]^2 Exp[χ[r]])/g[r]^2 φ[r] - m^2/g[r] φ[r] == 0;
eq2 = ϕ''[r] + (χ'[r]/2 + 2/r) ϕ'[r] - (2 qφ^2 φ[r]^2)/g[r]^2 ϕ[r] == 0;
eq3 = χ'[r] + r φ'[r]^2 + (r qφ^2 φ[r]^2 ϕ[r]^2 Exp[χ[r]])/g[r]^2 == 0;
eq4 = g'[r] + (1/r - χ'[r]/2) g[r] + (r ϕ'[r]^2 Exp[χ[r]])/4
- 3 r + (r m^2 φ[r]^2)/2 == 0;
eqns = {eq1, eq2, eq3, eq4, ϕ[ϵ] == ϕ0, ϕ'[ϵ] == ϵ,
φ[ϵ] == ϵ, φ'[ϵ] == ϵ, χ[ϵ] == χ0, g[ϵ] == ϵ};
{φ, ϕ, χ, g} = {φ, ϕ, χ, g} /. NDSolve[eqns, {φ, ϕ, χ, g}, {r, ϵ, 2}][[1]];
Plot[g[r], {r, ϵ, Sqrt[2]}]
Plot[ϕ[r], {r, ϵ, Sqrt[2]}]


• It works even without a minimal change. – István Zachar Dec 22 '15 at 20:25
• @István Zachar I know. Now we wait to see if he can uncover where in the editing process his non-working notebook turned into the working posted notebooks. – Bill Dec 22 '15 at 20:37
• Yes, eager to see the culprit! – István Zachar Dec 22 '15 at 21:47
• #Skylar15 There is something else hinting at, but not identifying, one or more underlying problems. If you change all the decimal constants to exact fractions so you can use higher precision and you add WorkingPrecision->128 and MaxSteps->10^6 (or variations of these) to NDSolve then the plots yank around and go off to infinity in various directions. Need to find the real underlying problem in your system. – Bill Dec 23 '15 at 18:16
• @Bill Sorry for the delayed response, but apparently the copy-paste is completely screwed up and 'phi' and 'varphi' are being swapped very inconsistently. Thanks for your effort though. I'm going to update the original code with correct 'phi', 'varphi'. I still get problem with that. – Skylar15 Dec 23 '15 at 18:24