I want to solve the following ODE:

```
y'(z)==-(y(z)^2-x(z)^2) chi/z^2
```

with the initial condition

```
y(z0) == x(z0)
```

where

```
x(z_) := -0.226679 E^(-0.991987 z) - 0.226679 E^(-0.991987 z) + 0.43999 E^(-0.965985 z);
chi = 5.5 10^12;
z0 = 20;
```

I know that the solution, i.e., y(z) should look like:

But instead of this, I got the next solution with Wolfram Mathematica (WM) using NDSolve:

My WM code is:

```
x(z_) := -0.226679 E^(-0.991987 z) - 0.226679 E^(-0.991987 z) +
0.43999 E^(-0.965985 z);
chi = 5.5 10^12;
z0 = 20;
solution =
NDSolve({ y'(z) == -(y(z)^2 - x(z)^2) chi/z^2, y(z0) == x(z0)},
y, {z, 10, 100},
Method -> {"EquationSimplification" -> "Residual"});
LogLogPlot({Legended(Evaluate(y(z) /. solution),
Placed(StyleForm("y", FontSize -> 12), {.7, .8})),
Legended(yxeq(z),
Placed(StyleForm("x", FontSize -> 12), {.7, .8}))}, {z, 100/7,
100}, PlotStyle -> {{Thickness(.004), Red}, {Thickness(.004),
Dashed, Purple}}, FrameLabel -> {"z", "y, x"},
AspectRatio -> 0.95,
PlotRange -> {{100/7, 100}, {0.2 10^-11, 4 10^-7}}, Frame -> True,
PlotStyle -> Thick,
FrameTicksStyle ->
Directive(FontSize -> 12, FontWeight -> Plain, FontColor -> Black),
LabelStyle ->
Directive(FontSize -> 12, FontWeight -> Plain, FontColor -> Black),
AspectRatio -> 0.95)
```

I do not understand why my solution is getting unstable. Could you help me with this issue?