Solution (source code)

= Solution
{c}

Put $\delta=\sqrt\epsilon$ and write the outer <asymptotic expansion> as $y_0+\delta y_{1/2}+\delta^2y_1+\cdots$. Equating powers in the <differential equation> gives
$$
y_0'-xy_0=0,\qquad y_{1/2}'-xy_{1/2}=0,\qquad y_1'-xy_1=y_0^{(3)}.
$$
The first two equations give $y_0=Ae^{x^2/2}$ and $y_{1/2}=Be^{x^2/2}$. Since $y_0^{(3)}=A(x^3+3x)e^{x^2/2}$, an <integrating factor> yields
$$
\boxed{y_0=Ae^{x^2/2},\quad y_{1/2}=Be^{x^2/2},\quad
 y_1=e^{x^2/2}\left[C+A\left(\frac{x^4}4+\frac{3x^2}2\right)\right].}
$$
Here $A,B,C$ are matching constants; boundary conditions have not yet been imposed.