Maybe replace the algebraic equation by a function definition. Reduce the AccuracyGoal
. I see no need for such a high setting, and it reduces the step size to zero. You could raise the PrecisionGoal
if a higher quality solution is needed.
Following integrates over the whole interval:
Block[{P1},
q[x_?NumericQ] := (830397. -
0.375 x^(87/296))/(4.51833*10^8 x^(34/37) - 1.02416*10^7 x -
1. x^(383/296));
P1[x_] := 2.372463129699238`*^-19 Sqrt[y2[x] + y4[x] + y1[x] + y3[x]];
sol1 = NDSolve[
SetPrecision[{y1'[x] == y1[x] (-3 P1[x] - q[x] + G'[x]/G[x]),
3 P1[x] y2[x] + y2'[x] == 0.95 y1[x] q[x],
4 P1[x] y3[x] + (y1[x] G'[x])/G[x] + y3'[x] == 0,
4 P1[x] y4[x] + y4'[x] == 0.05 y1[x] q[x],
(*P1[x]==2.372463129699238`*^-19 Sqrt[y2[x]+y4[x]+y1[x]+y3[x]],*)
G'[x] == (-8.722718060807123`*^73 +
1.80559635200202`*^-76 G[x]^4 y3[x])/G[x]^2,
3 y1[10^-30] == 10^51 \[Pi]^2, y2[10^-30] == 0,
y3[10^-30] == 3.3*10^56 \[Pi]^2, y4[10^-30] == 0,
G[10^-30] == 5.468912167778173`*^27}, 100], {y1, y2, y3, y4,
G}, {x, 10^-30, 10^5}, Method -> "Automatic", AccuracyGoal -> 20,
PrecisionGoal -> 5, WorkingPrecision -> 80]
]