Update: Incorporating comments
Originally I took for granted that the expected result was correct. A quick look at the IVP shows it starts on the nullcline $y(z)=x(z)$:
$$
y'(z)=-\frac{\chi\cdot\left(y(z)^2-x(z)^2\right)}{z^2}, \quad y(z_0)=x(z_0) \,.
$$
Since $y'(x)$ changes sign and $\chi$ is large, the quick reversal at $z=z_0$ should be expected.
The waviness to the right of $z=z_0=20$ in the OP is a numerical issue, but the issue is not stability. The default AccuracyGoal is around 8. Given an accuracy goal $a$, error estimates less than around $10^{-a}$ are accepted. When the magnitude of the solution is not much more or less than $10^{-a}$, then errors are accepted that are significantly large relative to the solution. I think of the accuracy goal as telling NDSolve the largest number that I would consider equivalent to zero. If it produces jittery noise that the user finds unacceptable, the user should increase AccuracyGoal. In this case, since the first run produced a solution that spends a lot of time around $10^{-11}$, one should expect to need an AccuracyGoal of around 19 or 20 to get single-precision results.
So two modifications are needed: Change the location of the initial condition to be to the left of the integration interval, and increase the AccuracyGoal.
However, being interested in the mathematical problem, I tried to see how close I could get to x[t]. We can do pretty well if we start at the maximum of x[t]. Since both y'[t] and x'[t] at this point, this is probably the best we can do; at all other initial conditions on the graph of x[t], the solution will cross the graph. The following does a good job until y[z] tries to change sign:
x[z_] = Rationalize[-0.226679 E^(-0.991987 z) -
0.226679 E^(-0.991987 z) + 0.43999 E^(-0.965985 z), 10^-6/2];
chi = 55/10*10^12;
z0 = z /. First@NSolve[x'[z] == 0 && 0 < z < 5, WorkingPrecision -> 32]
solution = NDSolve[
{y'[z] == -(y[z]^2 - x[z]^2) chi/z^2, y[z0] == x[z0]},
y, {z, 0, 100},
AccuracyGoal -> 20,
"ExtrapolationHandler" -> {Indeterminate &,
"WarningMessage" -> False}]
NDSolve::ndcf: Repeated convergence test failure at z == 1.1604817342549352`; unable to continue.
Changing variables $t=e^z$ and using the method "StiffnessSwitching", we can integrate further into y[z] < 0, but the integrator struggles when z < 0.7. The left edge of the plot below represents noise and not the true solution. A log plot is not possible because the solution changes sign, so I scaled the y coordinate by the varying factor Exp[z/2]. Exponential grid lines are shown with the scale indicated on top. The separation of y from x shows y[z] approaching $\approx 5.15 \times 10^{-12}$ as $z \rightarrow +\infty$.
z0 = z /. First@NSolve[x'[z] == 0 && 0 < z < 5, WorkingPrecision -> 32]
solution =
NDSolve[{Exp[-t] y'[t] == -(y[t]^2 - x[Exp[t]]^2) chi/Exp[t]^2,
y[Log@z0] == x[z0]}, y, {t, -Log@100, Log@100},
AccuracyGoal -> 20, PrecisionGoal -> 10,
Method -> "StiffnessSwitching",
"ExtrapolationHandler" -> {Indeterminate &,
"WarningMessage" -> False}, WorkingPrecision -> 30]
Show[
Plot[Evaluate[Exp[z/2] Join[
Table[10^y0/Exp[1/2], {y0, -12, -3}]]],
{z, 0, 43},
PlotStyle -> Join[
Table[Directive[AbsoluteThickness[0.5], Gray], {y0, -12, -3}]],
PlotRange -> {-0.015, 0.015}, Frame -> True,
FrameTicks -> {{Automatic, Automatic}, {Automatic,
Table[{1 - 2 (2 + y0) Log[10], 10.^y0}, {y0, -12, -3, 2}]}}],
Plot[Evaluate[Exp[z/2] Join[
{y[Log@z] /. solution, x[z]}]],
{z, 0, 42},
PlotStyle -> Join[
{AbsoluteThickness[3], AbsoluteThickness[1.5]}],
PlotRange -> All,
PlotLegends -> {y, x}]
]
It should be clear that the desired plot is somewhat misleading. On the one hand, the IC y[z0] == x[z0] always results in a solution that goes down to the left as in the OP. On the other, the plot makes it look like y[z] and x[z] are asymptotic to each other (or even that they agree) as z decreases toward zero.
Original answer
I think it's thought it was just numerical instability. [The analysis below the graph shows why the OP's plot should have a quick turn at the initial condition.] I had to change the y variable to $y = 10^{u/10}$, because StreamPlot is finicky and needs a plot range that is roughly square.
changeVar = y -> (10^(u[#]/10) &);
First@*Solve @@@ (
NestList[
D[#, z] /. First@Solve[First@#, u'[z]] &,
{y'[z] == -(y[z]^2 - x[z]^2) chi/z^2 /. changeVar, u'[z]}, 1]
) /. u[z] -> u // Flatten;
MapThread[Set[#[u_, z_], #2] &, {{up, upp}, {u'[z], u''[z]} /. %}];
?up
?upp
Show[
ContourPlot[
{up[u, z] == 0, upp[u, z] == 0},
{z, 15, 40}, {u, -115, -90},
ContourStyle -> {Directive[Dashed, Darker@Green],
Directive[Dashed, Darker@Magenta]},
FrameTicks -> {{Charting`ScaledTicks[{10 # &, #/10 &}],
Automatic}, {Automatic, Automatic}},
PlotLegends -> {u'[z] == 0, u''[z] == 0}
],
StreamPlot[
{1, up[u, z]},
{z, 15, 40}, {u, -115, -90},
StreamPoints -> {{{sepIC + 0.033, Red}, Automatic}}]
]
Moving to the left, if a solution steps across the green line u'[z] == 0, the solution heads down. Step across the purple line u''[z] == 0 and the solution becomes concave down and it will eventually step across the green line.
yxeq[z]as a function or you've meant to do something different than you do do. $\endgroup$AccuracyGoal -> 100toNDSolveimproves y for z>20 $\endgroup$AccuracyGoalat (64-bit) machine precision to be at least eight digits more than, say, the "average" value of the solution. $\endgroup$