Trying to find uniform approximation of a boundary value problem I have
$$ \varepsilon y'' + (1+x)^2 y' + y = 0 , x \in [0,1], y(0)=y(1)=1$$
I know I have a boundary layer at near $x=0$ since $(1+x)^2 > 0$. I want to find the uniform approximation of the solution up to order $\varepsilon^2$. for the outer region, using $y_{outer} \sim \sum y_n \varepsilon^n$, I obtained 
$$ y_{outer}(x) \sim e^{ \frac{1}{1+x} - \frac{1}{2} } + \varepsilon e^{   \frac{1}{1+x} - \frac{1}{2} } \left(  \frac{1}{2(1+x)^2} + \frac{1}{5(1+x)^3} -  \frac{3}{80} (1+x)^2 \right) $$
To find the inner solution, I introduce the scaling $X = x/\varepsilon$, and after some calculation, I obtained that 
$$ Y_{inner}(X) \sim A_0(1- e^{-X} ) - \varepsilon A_0 (  - (X^2+X+1) e^{-X} + X - B_0 e^{-X} + 1 + B_0 ) $$
now, I have to do the matching
For the leading order term matching, we must have 
$$ \lim_{x \to 0+} y_{outer} = \lim_{X \to \infty} Y_{inner} \implies e^{ \frac{1}{1+0} - \frac{1}{2} } = A_0( 1 - 0) \implies A_0 = e^{1/2} $$
however, Im stuck on fining the higher order matching. How is this done?
 A: Outer solution
For
\begin{align}
y_0&=e^{\frac1{1+x}-\frac12}\text{ and }y_1=uy_0
\end{align}
I get
\begin{align}
(1+x)^2y_1'+y_1&=-y_0''\\&=(1+x)^2y_0u'
\end{align}
so that, with $u(1)=0$,
\begin{align}
u'&=-\frac2{(1+x)^5}-\frac1{(1+x)^6}
\\
u&=\frac{1}{2(1+x)^4}+\frac1{5(1+x)^5}-\frac3{80}
\end{align}
This gives the first order approximation of the outer solution as
$$
y_{\rm outer}=e^{\frac1{1+x}-\frac12}+εe^{\frac1{1+x}-\frac12}\left(\frac{1}{2(1+x)^4}+\frac1{5(1+x)^5}-\frac3{80}\right)
$$
with
$$
y_{\rm outer}(0)=\sqrt{e}+ε\sqrt{e}\frac{53}{80}
$$

Inner solution
The inner equation in first perturbation order is $$Y''+(1+2εX)Y'+εY=0$$ where one demands $Y(0)=1-\sqrt{e}$, $Y(\infty)=0$ The zeroth order solution is $Y_0=C+De^{-X}$ which implies $C=0$ and $D=. 1-\sqrt{e}$. 
For the next perturbation order, extend the exponent to retain a limit at infinity. 
$$
Y=D·\exp(-aX-εbX^2),\;Y'=-(a+2εbX)Y,\;\\ Y''=(-2εb+(a+2εbX))^2)Y
$$
so that the ODE $Y''+(1+2εX)Y'+εY=0$ is satisfied to $O(ε^2)$ if 
\begin{align}
(-2εb+a^2)+(-a)+(ε)&=0\implies a=1+ε(2b-1)\text{ or }a=ε(1-2b)\\
4εab-2ε(a+b)X&=0\implies b=1+O(ε)\text{ or }b=0
\end{align}
The second variant is slow moving and can thus be subsumed by the outer solution. What remains is 
$$
Y(X)=D·\exp(-(1+ε)X-εX^2)
$$
and to second order error the coefficient is $D=1-\sqrt{e}(1+ε\frac{53}{80})$.
The full approximation is then
$$
y_{\rm approx}(x)=y_{\rm outer}(x)+Y(x/ε).
$$
