# Issue with NDSolve, Solving coupled differential equations, Error test failure at x = some value

yEQ[x_] :=
45/(4*3.14^4) gChi/gs x^2 BesselK[2, x] /. {gChi -> 2,
gs -> 106.7, MP -> 2.44*10^18(*GeV*)}
A = (2 Sqrt[2 ]*3.14)/(3 Sqrt) gs^(1/2) MP /. {gChi -> 2,
gs -> 106.7, MP -> 2.44*10^18(*GeV*)};
init := 0.01

SolveBEQNN[mChi_, SigmaV_] =
NDSolve[{y'[x] == (-A mChi)/x^2 SigmaV*(y[x]^2 - (yEQ[x])^2),
y[init] == (yEQ[init])}, y, {x, init, 100000},
InterpolationOrder -> All, MaxSteps -> 10000000,
MaxStepFraction -> 1/100000, AccuracyGoal -> \[Infinity]];

PSPI = y[2*50000] /. SolveBEQNN[100, 1.7*10^-9.9]


The above code gives correct result, however when I change the magnitude of the variables e.g.

 PSPI = y[2*50000] /. SolveBEQNN[100, 1.7*10^-8]
NDSolve::nderr: Error test failure at x == 0.5118756600640246; unable to continue.
InterpolatingFunction::dmval: Input value {100000} lies outside the range of data in the
interpolating function. Extrapolation will be used.


it is giving the following output with error messages.

Out: {-1.14979*10^15}


which is not correct. How to avoid this and find the correct solution.

• Your code in any case gives the wrong result. – Alex Trounev Mar 17 '20 at 12:05
• Can you specify the reason? – Abhijit Saha Mar 17 '20 at 14:36
• I just ran your code and got three messages: NDSolve::ndnum: Encountered non-numerical value for a derivative at x == 0.01. NDSolve::nderr: Error test failure at x == 0.03051719479497417; unable to continue. InterpolatingFunction::dmval: Input value {100000} lies outside the range of data in the interpolating function. Extrapolation will be used. – Alex Trounev Mar 17 '20 at 14:42
• Please check with PSPI = y[2*50000] /. SolveBEQNN[100, 1.7*10^-11], It will provide results without showing the last two error message. First one is coming because I am using variables in NDSolve. – Abhijit Saha Mar 17 '20 at 14:46
• Do you want to check your code or solve a problem? – Alex Trounev Mar 17 '20 at 17:06

I would definitely have opted for ParametricNDSolve for equations with unknown parameters; NDSolve, I use it for the equations which with deterministic variables. Also, Set delayed can often cause process to fail. Another change I made is using StiffnessSwitching method. This is what I did:

    yEQ[x_] =  45/(4*3.14^4) gChi/gs x^2 BesselK[2, x] /. {gChi -> 2, gs -> 106.7, MP -> 2.44*10^18(*GeV*)}
A = (2 Sqrt*3.14)/(3 Sqrt) gs^(1/2) MP /. {gChi -> 2, gs -> 106.7, MP -> 2.44*10^18(*GeV*)};
init = 0.01


Now

    SolveBEQNN = ParametricNDSolve[{y'[x] == (-A mChi)/x^2 SigmaV*(y[x]^2 - (yEQ[x])^2),
y[init] == (yEQ[init])}, y, {x, init, 100000}, {mChi, SigmaV},  Method -> "StiffnessSwitching"]


gives

   y[100, 1.7*10^-8][2*50000] /. SolveBEQNN
(*7.18485*10^-12*)

• Thanks a lot, this works. – Abhijit Saha Mar 18 '20 at 5:47