Solution (source code)

= Solution

Write $E=\operatorname{erf}z$, $g=e^{-z^2}$ and $b=2/\sqrt\pi$. The first correction in the <inner expansion> satisfies $\mathcal L_0Y_1=2zE$. Direct differentiation verifies that $zE+bg$ is a <particular solution>. Adding the constant and <error function> <homogeneous solutions> gives $Y_1=zE+bg+A+B E$.

The left <boundary condition> forces $A=-b$. Re-expanding the <outer expansion> at $x=sz$ gives
$$
y_{\mathrm{out}}(sz)=1+sz+s^2\left[\frac{z^2}{2}-\frac12\log z-\frac12\log s\right]+\cdots.
$$
Therefore $Y_1\sim z$, with no constant term in its large-$z$ overlap. Since $E\to1$ and $g\to0$, this sets $B=b$. Thus the required half-order correction is
$$
\boxed{Y_1(z)=z\operatorname{erf}z+\frac2{\sqrt\pi}\left[e^{-z^2}+\operatorname{erf}z-1\right].}
$$
It vanishes at zero and has the correct linear overlap. Keeping this $O(\sqrt\varepsilon)$ term is essential even though the outer solution has no correction of that order.