Roughly speaking, an AccuracyGoal of
$g$ means absolute errors less than
$10^{-g}$ are accepted/acceptable. When the magnitude of the solution gets less than this threshold, then errors greater in magnitude may cause the solution to change sign. I think of AccuracyGoal as telling NDSolve the threshold below which a solution may safely be treated as zero plus or minus an error. In this case, you don't want that. Theoretically raising AccuracyGoal so that
$10^{-g}$ is much less than the magnitude of the solution is a fix. However, the corrector step in LSODA sometimes/often will exhibit a convergence problem, which it did for me. I could get that to go away by keeping PrecisionGoal modestly low.
So here's the first solution:
sol12a = NDSolve[{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),
y1[10] == 10^-7, y2[10] == 0}, {y1, y2}, {x, 10, 10^12},
Method -> "Automatic", PrecisionGoal -> 6,
WorkingPrecision -> 100, MaxSteps -> \[Infinity]];
Now, the first equation is self-contained subsystem with an exact solution. So here's a second method that avoid the numerical issue entirely:
y1sol = DSolve[{y1'[x] + 4 y1[x]/x == 0, y1[10] == 10^-7}, y1, x];
y2sol =
NDSolve[{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),
y2[10] == 0} /. First@y1sol,
{y2}, {x, 10, 10^12}, Method -> "Automatic",
WorkingPrecision -> 30, MaxSteps -> \[Infinity]];
sol12b = Join[y1sol, y2sol, 2];
I usually do something like Daniel's approach, but in this case I'd scale the variable y1, which is the misbehaving one. Transforming
$x=e^t$ might also be convenient — or not. It often is when a variable ranges from
$10$ to
$10^{12}$. Here's a third solution:
eqns0 = {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)};
eqns1 = eqns0 /. y1 -> Exp@*u1 // Simplify[#, u1[x] > 0] &;
sol12c = NDSolve[Join[eqns1, {u1[10] == Log[10^-7], y2[10] == 0}],
{u1, y2}, {x, 10, 1*^12}] /. (u1 -> usol_) :> y1 -> Exp@*usol;
