Solving a stiff nonlinear ODE system

The system I am trying to solve is simple, but looks pretty stiff and I have unsuccessfully tried to solve it with StiffnessSwitching. It is the following one:

x0=.001;
y0=10^64;
const=3*10^9;
h[a_] := Sqrt[y[a] + x[a]];
Needs["DifferentialEquationsNDSolveProblems"];
Needs["DifferentialEquationsNDSolveUtilities"];
system = {
a*y'[a] + 3*y[a] == -const*y[a]/h[a],
a*x'[a] + 4*x[a] == const*y[a]/h[a],
x == x0,
y == y0
};

system2 = ReplacePart[system /. {x -> (F[#]^2 &), y -> (G[#]^2 &)}, {{3} -> F == Sqrt[x0], {4} -> G == Sqrt[y0]}];
{nds1, nds2} = NDSolveValue[system2, {F, G}, {a, 1, 10^7}, Method -> "StiffnessSwitching", MaxSteps -> 10^6, WorkingPrecision -> MachinePrecision];

where system2 has been defined to ensure the positiveness of the functions. But here I receive the stiffness message:

NDSolveValue::ndsz: At a == 1.`, step size is effectively zero; singularity or stiff system suspected.

I have also tried with ExplicitRungeKutta, but nothing happened. It is clear that the initial condition y0 is extremely high, but I need it. Trying with a lower y0=10^16 it works, but only if I do not include const.

Any hint?

• Evaluating system2[[1;;2]] at a == 1 yields {F' -> 4.74342*10^42, G' -> -1.5*10^32}. Mathematica probably concludes that these enormous values indicate stiffness. Try rescaling your variables to obtain smaller coefficients. By the way, you do not need the code containing Needs. Sep 2 '19 at 19:14
• Which coefficients are you referring to?
– Lele
Sep 4 '19 at 14:43
• y0 in particular. Sep 4 '19 at 14:48
• Yes, with lower y0 it works, as I said. But, physically, I need that high value!
– Lele
Sep 4 '19 at 15:27
• Try rescaling your variables so that these numbers are not so large. Sep 4 '19 at 15:47

The computation can be performed as follows. First, solve for {x, y} instead of {F, G}, because the ODEs are simpler, and then obtain {F, G} by taking the square roots of {x, y}. Even then,

NDSolveValue[system, {x, y}, {a, 1, 10^7}]

immediately fails, claiming that the "step size is effectively zero; singularity or stiff system suspected". Perhaps, NDSolve concludes this, because the second equation in system is approximately a*x'[a] + 4*x[a] == 3 10^41. Rescaling dependent variables to eliminate such large constants. i.e., y by y0 and x by const Sqrt[y0], probably would allow the computation to proceed. In addition, changing the independent variable from a to Log[a] reduces the enormous domain of integration to a more manageable size. As it turns out, this last transformation alone is sufficient to obtain the solution. With b == Log[a], system becomes

system1 = {y'[b] + 3*y[b] == -const*y[b]/h[b], x'[b] + 4*x[b] == const*y[b]/h[b],
x == x0, y == y0};

which is, incidentally, autonomous. Then,

NDSolveValue[system1, {x[b], y[b]}, {b, 0, Log[10^7]}];
LogPlot[%, {b, 0, Log[10^7]}, ImageSize -> Large,
AxesLabel -> {b, "{x,y}"}, LabelStyle -> {Bold, Black, 15}]
Plot[First@%%, {b, 0, 2}, PlotRange -> All, ImageSize -> Large,
AxesLabel -> {b, x}, LabelStyle -> {Bold, Black, 15}]  The second plot, a blowup of x near the origin, shows that it grows rapidly but linearly there. Why transforming only the independent variable is sufficient to obtain a solution is unclear. Perhaps, NDSolve uses a different computational method for autonomous ODEs.

• Great! Thank you! It was much easier than I had thought... last question: how do I use the solutions from NDSolveValue? I tried with {xsol,ysol}=NDSolveValue[system1, {[b], y[b]}, {b, 0, Log[10^7]}];, but when I evaluate one of them, for example, as xsol it returns only InterpolatingFunction[...][b]. I would like to use those solutions to define other functions and so on... thank you!
– Lele
Sep 5 '19 at 10:27
• Try xsol /. b -> 10. Alternatively, start with {xsol, ysol} = NDSolveValue[system1, {x, y}, {b, 0, Log[10^7]}]; (Look at the output from NDSolveValue using {x[b], y[b]} vs {x, y} to see the difference.) Sep 5 '19 at 11:50