# Stiffness problem in an NDSolve system. StiffnessSwitching does not help?

I'm trying to solve an ODE system with the NDSolve method. This my ODE system, with the BC and the functions V and dV defined: s = NDSolve[{F'[
t] == -(2 /Sqrt) (F[t]^2/2 + a[t]^2*  V[B[ t]])^(1/2)*F[t] -
a[t]^2*dV[B[t]],
a'[t] == a[t]/Sqrt (F[t]^2/2 + a[t]^2*  V[B[ t]])^(1/2),
B'[t] == F[t], F[tmin] == Subscript[v, 0], a[tmin] == a0 ,
B[tmin] == Subscript[B, 0] }, {F, a, B}, {t, tmin, tmax}, Method->"ExplicitRungeKutta"]


Where the solution functions are:

F(t), a(t), B(t)


I need the solutions for a big interval, like:

t = [-10^(10), -10^(-10)]


So:

tmin = -10^(10)
tmax = -10^(-10)


But the system seems to be stiff. I'm trying to replace the "ExplicitRungeKutta" method for "StiffnessSwitching method", replacing the method with:

Method -> {"StiffnessSwitching", Method -> {"ExplicitRungeKutta", Automatic}}


But the system seems to remain stiff.

The complete code is:

Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
Needs["NumericalCalculus"];

V[B_] := 3;
dV[B_] :=  V'[B];

B0 = 1;
F0 = 0;
a0 = 10^(-10);

tmin = -10^(10);
tmax = -10^(-10);

s = NDSolve[{F'[
t] == -(2 /Sqrt) (F[t]^2/2 + a[t]^2*  V[B[ t]])^(1/2)*F[t] -
a[t]^2*dV[B[t]],
a'[t] == a[t]/Sqrt (F[t]^2/2 + a[t]^2*  V[B[ t]])^(1/2),
B'[t] == F[t], F[tmin] == F0, a[tmin] == a0 , B[tmin] == B0 }, {F,
a, B}, {t, tmin, tmax},  Method->"ExplicitRungeKutta"]


Question: Is there any way to obtain the solutions of these equations without stiffness? What method should I use?

I'd appreciate much all the answers.

• What if you just use Method -> "StiffnessSwitching" without invoking RK at all? – J. M.'s discontentment Jun 21 '16 at 23:03
• (1) What's v0? (2) What if it's a singularity and not stiffness? – Michael E2 Jun 22 '16 at 1:02
• Are tmin & tmax switched, or is tmin supposed to be greater than tmax? – Michael E2 Jun 22 '16 at 3:09
• V0 was F0, sorry for the mistake. And both variables tmin and tmax were negative (another mistake writting the code here). – Guillermo Martínez Somonte Jun 22 '16 at 23:57
• And the stiffnessSwitching method alone does not solve the (aparent) stiffness problem. – Guillermo Martínez Somonte Jun 22 '16 at 23:59

NDSolveIterate::ndsz: At t == -1.*10^10, step size is effectively zero; singularity or stiff system suspected.

The "step size is effectively zero" means that in floating-point arithmetic t + dt is equal to t for the computed time step dt. If this is the problem, increasing working precision might help.

de = {F'[t] == -(2/Sqrt) (F[t]^2/2 + a[t]^2*V[B[t]])^(1/2)*F[t] -
a[t]^2*dV[B[t]],
a'[t] == a[t]/Sqrt (F[t]^2/2 + a[t]^2*V[B[t]])^(1/2),
B'[t] == F[t]};
ics = {F[tmin] == F0, a[tmin] == a0, B[tmin] == B0};
s =
NDSolve[{de, ics}, {F, a, B}, {t, tmin, tmax},
PrecisionGoal -> 8, AccuracyGoal -> 8, WorkingPrecision -> 20,
Method -> "StiffnessSwitching"] Note: The first of the OP's differential equations has

F'[t] == <..stuff..> * F[t]


because dV[B] is zero and the second terms goes away. Since F[t] starts out at F[tmin] == 0, that means F[t] is a constant 0 (zero). Thus the third differential equation B'[t] == F[t] and its initial condition imply B[t] is a constant 1 (one).