ODE: Asymptotic inner expansions I am given the following ODE $$ \epsilon \frac{d^2y}{dx^2}+x\frac{dy}{dx}+y=0 \; , \qquad 1<x<2$$
with the BVPs $$y(1)=0 \qquad , \qquad y(2)=1$$
and where $0<\epsilon \ll1$.
I am to use the method of matching inner and outer expansions to solve this.
By letting $y(x,\epsilon) \sim y_0(x)+y_1(x)\epsilon + y_2(x)\epsilon ^2+\cdots$, I have found the leading-order term of the outer expansion to be $2/x$ (and this satisfies only $y(2)=1$).
In computing the inner expansion, I let $x=\epsilon ^a z$ so that we now have $$\epsilon ^{1-2a}\frac{d^2Y}{dz^2}+z\frac{dY}{dz}+Y=0$$
where $Y(z,\epsilon)=y(x,\epsilon)$.
In order to match the orders, we need to choose $a$ such that $1-2a=0$, thus $a=\frac 12$. Plugging this in, we have $$\frac{d^2Y}{dz^2}+z\frac{dY}{dz}+Y=0$$
Letting $Y(z,\epsilon) \sim Y_0(z)+\cdots$, the equation at leading order is simply $$\frac{d^2Y_0}{dz^2}+z\frac{dY_0}{dz}+Y_0=0$$ and solving this gives $$Y_0=e^{-x^2/2}\Bigl(A+B \; \text{erfi} \bigl(\frac{x}{\sqrt 2} \bigr) \Bigr)$$
which doesn't seem to be right. Have I done something wrong? Or, ignoring what I have done, what is the right way to approach this?
 A: For the fast movements of the solution you in general look for local coordinates $x=x_0+ϵ^aX$ which gives for $Y(X)=y(x_0+ϵ^aX)$
$$
ϵ^{1-2a}Y''+(x_0+ϵ^aX)ϵ^{-a}Y'+Y=0\iff ϵ^{1-a}Y''+(x_0+ϵ^aX)Y'+ϵ^{a}Y=0
$$
which gives $a=1$ for the width exponent of the inner solution (as the exceptional point $x_0=0$ is not inside the integration interval). As $x_0>0$ the solution to the equation of the leading terms $Y''+x_0Y'=0$ is $Y(X)=Ce^{-x_0X}+D$. As it is not bounded in direction $X\to-\infty$ the boundary layer has to occur at the left interval end $x_0=1$ with first order approximation of the solution as
$$
y(x)=\frac2x -2e^{-(x-1)/ϵ}
$$

$ϵ=0.01$, blue = numerical solution, green = boundary layer approximation
A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$

Lets $\ds{\,\mrm{y}\pars{x} \equiv
\exp\pars{\mrm{S}\pars{x,\delta} \over \delta}\,,\quad
\delta \equiv \root{\epsilon}.\quad}$

$$
\mbox{Then,}\quad\left\{\begin{array}{rcl}
\ds{\mrm{y}'\pars{x}} & \ds{=} &
\ds{\mrm{y}\pars{x}{\mrm{S}'\pars{x,\delta} \over \delta}}
\\
\ds{\mrm{y}''\pars{x}} & \ds{=} &
\ds{\mrm{y}\pars{x}{\bracks{\mrm{S}'\pars{x,\delta}}^{\,2} \over \delta^{2}} +
\mrm{y}\pars{x}{\mrm{S}''\pars{x,\delta} \over \delta}}
\end{array}\right.
$$
and
\begin{align}
&\delta^{2}\braces{\mrm{y}\pars{x}{\bracks{\mrm{S}'\pars{x,\delta}}^{\,2} \over \delta^{2}} +
\mrm{y}\pars{x}{\mrm{S}''\pars{x,\delta} \over \delta}} +
x\bracks{\mrm{y}\pars{x}{\mrm{S}'\pars{x,\delta} \over \delta}} + \mrm{y}\pars{x} = 0
\\[5mm] &
\bbx{\bracks{\mrm{S}'\pars{x,\delta}}^{\,2}\,\delta +
\mrm{S}''\pars{x,\delta}\delta^{2} +
x\,\mrm{S}'\pars{x,\delta} + \delta = 0}
\end{align}

Now, you can expand $\ds{\,\mrm{S}\pars{x,\delta}}$ in powers of $\ds{\delta}$. 

Namely,
$$
\mrm{S}\pars{x,\delta} = \mrm{S}_{0}\pars{x} + \mrm{S}_{1}\pars{x}\delta  + \mrm{S}_{2}\pars{x}\delta^{2} + \cdots  
$$
$$
\left\{\begin{array}{rcl}
\ds{\mrm{S}_{0}'\pars{x}} & \ds{=} & \ds{0}
\\[1mm]
\ds{x\,\mrm{S}_{1}'\pars{x} + 1} & \ds{=} & \ds{0}
\\[1mm]
\ds{\mrm{S}_{2}'\pars{x}} & \ds{=} & \ds{0}
\\[1mm]
\ds{\mrm{S}_{1}'^{2}\pars{x} + \mrm{S}_{1}''\pars{x} + x\,\mrm{S}_{3}'\pars{x}} & \ds{=} & \ds{0}
\\[1mm]
\ds{\vdots} & \ds{\vdots} & \ds{\vdots}
\end{array}\right.
$$

This approach is the
  WKB Approximation.

A: Doing the WKB approximation right
As Felix Marin cited, one starts by representing the solution as $y(x)=\exp\left[\frac1δS(x;δ)\right]$ where later on $S(x;δ)$ is developed as power series in $δ$, $S(x;δ)=S_0(x)+S_1(x)δ+S_2(x)δ^2+\dots$. Inserting the derivatives $y'=\frac1δyS'$ and $y''=\frac1δyS''+\frac1{δ^2}yS'^2$ into the differential equation gives
$$
0=ϵ\left(\frac1δyS''+\frac1{δ^2}yS'^2\right)+\frac1δxyS'+y
=\frac1δy\left(\fracϵδS'^2+ϵS''+xS' +δ\right).
$$ 
In the last form one sees that the dominant terms in the last factor are 
$$
\fracϵδS'^{\,2}+xS'.
$$ 
To get meaningful results from the approximation method these terms have to be balanced in magnitude. This happens for $δ=ϵ$ with the resulting equations for the perturbation series 
\begin{align}
&&(S'^2+xS')+δ(S''+1)&=0\\[1em]\hline
&δ^0:& S_0' (S_0' +x)&=0
\\
&δ^1:& 2S_0' S_1'  +xS_1'  +S_0''+1&=0
\\
&δ^2:& S_1'^2+2S_0'S_2'+xS_2'+S_1''   &=0
\\
&&&\vdots
\end{align} 
The first equation has two different solutions for $S_0'$ that give rise to two basis solutions for $y$. With $S_0'=-x$ the following equations result in
$S_1' =S_2' =0,...$. The second solution has $S_0' =0$, thus no fast term, and leads to $S_1' =-\frac1x$, $S_2' =-\frac2{x^3},...$.
Integration of the terms leads to the assembled general solution approximation
$$
y(x)=A\exp\left(-\frac{x^2-1}{2ϵ}+O(ϵ^2)\right)+B\frac1x\exp\left(\frac{ϵ}{x^2}-ϵ+O(ϵ^2)\right),
$$ where integration constants were chosen to give $0=y(1)=A+B$. The second boundary condition then gives
$$
1=y(2)=A\exp\left(-\frac{3}{2ϵ}\right)+\frac{B}2\exp\left(-\frac34ϵ\right)
\implies
-A=B=2\exp\left(\frac34ϵ\right)
$$
ignoring the very small term $\exp(-\frac{3}{2ϵ})$.

Comparing against the numerical solution confirms a good fit even for $ϵ$ as large as $0.1$, while moving away visibly for larger $ϵ$.

$ϵ=0.1$, blue = numerical solution, red = $O(ϵ^2)$ WKB approximation


$ϵ=0.2$, blue = numerical solution, red = $O(ϵ^2)$ WKB approximation
