ODE with nested boundary layers Problem:
Consider the equation
$$\varepsilon^3 \frac{d^2y}{dx^2} + 2x^3 \frac{dy}{dx} - 4\varepsilon y = 2x^3 \qquad \qquad y(0) = a \;, \; y(1)=b$$
in the limit as $\varepsilon \rightarrow 0^+$, where $0<a<b-1$.
Assuming that there are nested boundary layers at $x=0$, determine the thickness of the boundary layers and the leading-order additive composite solution.

Question:
My attempt is shown below. Basically, I don't understand how you are meant to match the different solutions in the various domains.
Would be grateful if someone could explain this to me, or give me some hints.

Attempt:
Firstly, the leading term $y_0$ of the outer solution satisfies
$$2x^3 \frac{dy_0}{dx} =2x^3 \qquad \qquad y_0(1)=b$$
This is easily solved to give $\color{red}{y_0(x) = x+b-1}$.
Now suppose we scale the equation with $x=\varepsilon ^\alpha X$ where $\alpha>0$ and $X = \mathcal O(1)$. The equation becomes
$$\varepsilon^{3-2\alpha} \frac{d^2y}{dX^2} + 2\varepsilon^{2\alpha}X^3 \frac{dy}{dX} - 4\varepsilon y = 2\varepsilon^{3\alpha}X^3$$
The possible leading order balances come from $\alpha = 1$ and $\alpha = 1/2$.
When $\alpha = 1/2$ (i.e. $x = \varepsilon^{1/2} X$), the leading order term $Y_0$ in this layer should satisfy
$$2X^3 \frac{dY_0}{dX} - 4Y_0 = 0 \qquad \qquad Y_0(X=0) = a$$
The general solution is $\color{red}{Y_0(X) = A\exp (-1/X^2)}$, which cannot satisfy the boundary condition, since $a>0$. So I suppose this is why we need another boundary layer.
Scaling instead with $x = \varepsilon \tilde X$, the leading order solution $\tilde Y_0$ in this layer satisfies
$$\frac{d^2\tilde Y_0}{d\tilde X^2} - 4\tilde Y_0=0\qquad \qquad \tilde Y_0(\tilde X = 0) = a$$
The solution is $\color{red}{\tilde Y_0(\tilde X) = \tilde A\sinh(2\tilde X)} \color{green}{+a\cosh(2\tilde X)}$.
$\color{blue}{\text{Now the problem is, how do I match these three solutions?}}$
From what I have learnt, I need to do something like
$$\lim_{x\rightarrow 0^+} y_0(x) = \lim_{X\rightarrow +\infty}Y_0(X)$$
to obtain $\color{blue}{A= b-1?}$ But if I similarly try to do
$$\lim_{X\rightarrow 0^+}Y_0(X) = \lim_{\tilde X\rightarrow +\infty}\tilde Y_0(\tilde X)$$
it doesn't work, because $\tilde Y_0$ is not bounded??
 A: You have three solutions to work with, the outer solution
$$ y_0(x) = x+B $$
the "wide" inner solution
$$ Y_0(X) = C_1\exp\left(-\frac{1}{X^2}\right),\quad X=\frac{x}{\sqrt{\epsilon}} $$
and the "narrow" interior solution
$$ \tilde Y_0(\tilde X) = D_1\exp\left(-2\tilde X\right)+D_2\exp\left(2\tilde X\right),\quad\tilde X=\frac{x}{\epsilon}. $$
Immediately you can say $B=b-1$, $D_2=0$ and $D_1=a$. 
There's no need to use $\cosh$ and $\sinh$, they're just combinations of $\exp(2X)$ and $\exp(-2X)$ anyway. The whole reason you set $D_2=0$ is to ensure that $\tilde Y_0(\tilde X)$ is bounded as $\tilde X\to\infty$, but $\cosh$ and $\sinh$ are both unbounded.
You have two places where you need to do asymptotic matching (to find constant or verify things match), which you identified correctly,
$$ \lim_{X\to\infty}Y_0(X)=\lim_{x\to0^+}y_0(x) \Rightarrow C_1=b-1,\qquad(\star)$$
and 
$$ \lim_{X\to0^+}Y_0(X)=\lim_{\tilde X\to\infty}\tilde Y_0(\tilde X) \Rightarrow 0=0 \qquad(\dagger),$$
so the inner solutions are exponentially small in the transition region.
You can combine all this to get a uniform approximation,
$$y(x) = \underbrace{\left[x+b-1\right]}_\text{outer}+\underbrace{\left[(b-1)\exp\left(-\frac{\epsilon}{x^2}\right)\right]}_\text{wide inner}+\underbrace{\left[a\exp\left(-\frac{2X}{\epsilon}\right)\right]}_\text{narrow inner}-\underbrace{\left[b-1\right]}_\text{matching constant $(\star)$}-\underbrace{\left[0\right]}_\text{matching constant $(\dagger)$}+O(\epsilon^{1/2}),$$
that is,
$$y(x) = x+(b-1)\exp\left(-\frac{\epsilon}{x^2}\right)+a\exp\left(-\frac{2\color{red}{x}}{\epsilon}\right)+O(\epsilon^{1/2}),$$
A: Illustration of the solution above: A numerical solution of the boundary value problem is plotted against the approximation to show that indeed the 3 parts do fit together that way for sufficiently small $ε$, visually for $ε<0.05$.

