Here’s a way which shows off the "Splitting" and "LocallyExact" methods of NDSolve[]. I rarely use "Splitting" because it requires that one can split the vector field of the ODE f=f_1+f_2 in a nice way as well as, in practice, that the ordinary integration methods fail. Consequently I have to look things up to remind myself how it is implemented. It seems easier to me to use DSolve[] as I did in my first reply. But I wanted to remind myself how splitting works. Since "LocallyExact" uses DSolve[] under the hood to advance time integration, I was curious to see if I could get NDSolve[] to do it all for me.
"Splitting" is explained here. I followed an example there. One requirement is that the ODE be autonomous. I used Daniel’s substitution x=e^t in the form x == z[t], z'[t] == z[t], z[Log[10]] == 10, which will convert John Wick’s nonautonomous system of two DEs into an autonomous system of three DEs. One can then split the vector field into a sum of the form
\cases{y_1' = -4 y_1 &\cr y_2' = 0 & \cr z' = z &\cr} + \quad\cases{y_1' = 0 &\cr y_2' = *** & \cr z' = 0 &\cr}
The first can be integrated with "LocallyExact". For the second, we’ll use one-step numerical method. "Splitting" seems to use a fixed step size, which you can set with StartingStepSize, or you can let NDSolve[] estimate it for you (presumably based on the order and error estimates at the beginning of the integration). If you want NDSolve[] to check and reset the step size periodically, I think adding an event like WhenEvent[Mod[t, 1] == 0, {t; "RestartIntegration"}] works. The no-op t seems necessary to trick NDSolve[] into recalibrating.
P1 = 10^7;
P2 = 10^15;
y3 = 10^-6*(E^-x)*x^(-3/2);
ode = {y1'[x] + 4 y1[x]/x == 0,
y2'[x] + (3/x) (((\[Sqrt](y1[x] + y2[x]))/((\[Sqrt](y1[x] + y2[x]) -
P1 ((y2[x] - y3)/y1[x])) -
P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[x] == (-4 P2 (y2[x]^2 - y3^2) -
4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) -
P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x)};
ode = Join[
ode /. {f_'[x] :> f'[t]/z[t], f_[x] :> f[t], x -> z[t]} //
Simplify[# // Map[MultiplySides[#, z[t]] &], z[t] > 0] &,
{z'[t] == z[t]}];
ode = ode /. {f_'[t] :> f'[t], f_[t] :> f[t], t -> z[t]};
(*Print[ode/.v_[t]:>v];*)
With[{sol = First@Solve[ode, {y1'[t], y2'[t], z'[t]}]},
f1 = {y1'[t] == (y1'[t] /. sol), y2'[t] == 0,
z'[t] == (z'[t] /. sol)};
f2 = {y1'[t] == 0, y2'[t] == (y2'[t] /. sol), z'[t] == 0};
];
(*Print[{f1,f2}];*)
sol12 = NDSolve[
{ode, y1[Log[10]] == 10^-7, y2[Log[10]] == 0, z[Log[10]] == 10}
, {y1, y2, z}, {t, Log[10], Log[10^12]}
Method -> {"Splitting", "DifferenceOrder" -> 6,
"Equations" -> {f1, f2, f1}
, "Method" -> {
"LocallyExact" (* f1 *)
, {"ImplicitRungeKutta" (* f2 *)
, "Coefficients" -> "ImplicitRungeKuttaGaussCoefficients"
, "DifferenceOrder" -> 6
, "ImplicitSolver" -> {FixedPoint,
AccuracyGoal -> MachinePrecision, PrecisionGoal -> MachinePrecision,
"IterationSafetyFactor" -> 1}}
, "LocallyExact"}} (* f1 *)
, WorkingPrecision -> MachinePrecision
(*,InterpolationOrder->All*)
]
p1 = LogLogPlot @@ {{y1[t]/(y1[t] + y2[t]), y2[t]/(y1[t] + y2[t])} /.
t -> Log[x] /. sol12, Flatten@{x, Exp[y1["Domain"] /. sol12]},
PlotRange -> Full}
