Asymptotic approximation for $\epsilon^2 y'' = axy$ 
I would like to find the first two terms of the asymptotic approximation for $\epsilon^2 y'' = axy$ on $0\leq x <
 \infty$ where $y(\infty) = 0$ and $a>0$.

Here is my work so far, with the standard WKB method, we plug in 
$$y(x) = \exp\left[\frac{1}{\delta} \sum_{n=0}^\infty \delta^n S_n(x)\right]$$ and divide by it; our ODE simplifies to 
$$\epsilon^2\left[\frac{1}{\delta} S_0'' + S_1 + \cdots + \frac{1}{\delta^2} S_0'^2 + 2 \frac{1}{\delta} S_0'S_1' + \cdots \right] -ax = 0 \quad \quad (\star)$$
Our largest term $\frac{\epsilon^2}{\delta^2} S_0'^2 $must match with $-ax$, so we make $\delta = \epsilon$, and get 
$$S_0'^2 \sim ax$$
thus 
$$S_0 \sim \pm \sqrt{a} \frac{2}{3} x^{3/2}.$$
And by matching $O(\epsilon)$ terms in $(\star)$, we get 
$$2S_0'S_1' + S_0'' = 0$$
and we get $S_1 \sim - \frac{1}{4} \log(ax)$.
So our approximation 
$$y (x)\sim C_1 (ax)^{-1/4} \exp \left[\frac{1}{\epsilon} \sqrt{a} \frac{2}{3} x^{3/2}\right] + C_2 (ax)^{-1/4} \exp \left[-\frac{1}{\epsilon} \sqrt{a} \frac{2}{3} x^{3/2}\right]$$
And from our condition $y(\infty) = 0$, we get $C_1 = 0$. 
But the above approximation does not make sense at $x=0$, how can we take care of this problem? 
I know boundary layer theory but often it has problems when the coefficient of $y'$ is constant zero. I also worked out the WKB method for a boundary layer problem here, I am not sure how would I apply to this case.
Edit: If I add a boundary layer at $x=0$, we let $\xi = x/\delta$, then the ODE in terms of delta will become 
$$\frac{\epsilon^2}{\delta^2} (y_0'' + \epsilon y_1'' + \cdots ) = a\xi\delta(y_0 + \epsilon y_1, + \cdots)$$
which gives 
$$\frac{\epsilon^2}{\delta^2} y_0'' + \frac{\epsilon^3}{\delta^2} y_1'' + \cdots  - \delta a\xi y_0 + \cdots = 0$$
If I let $\delta = \epsilon$, I will get $y_0 = a\xi + b$, how would I match this with $y(x)$ from WKB method
 A: Your asymptotic approximation is correct :
$$y (x)\sim C_2 (ax)^{-1/4} \exp \left[-\frac{1}{\epsilon} \sqrt{a} \frac{2}{3} x^{3/2}\right] \tag 1$$
according to the boundary condition $y(\infty)=0$.
Note : A second order ODE requires at least two boundary conditions to expect having a unique solution. Since the wording of the question specifies only the condition $y(\infty)=0$, the solution remains undetermined : The coefficient $C_2$ cannot be computed.
Another key point is that an asymptotic formula is valid only for $x\to\infty$, not for $x\to 0$.So, don't expect Eq.$(1)$ makes sens at $x=0$. It is normal that the asymptotic approximation be far to be correct for small values of $x$.
If you want an approximate around $x=0$, expand $y(x)$ in power series. 
$$y(x)\simeq c_0+c_1x+c_2x^2+...$$
Of course, this formula is not valid for large $x$, so the boundary condition at infinity cannot be fulfilled. Moreover, two boundary conditions should be required close to $x=0$ in order to determine the coefficients $c_0$ , $c_1$, $c_2$ , …
In fact, as Claude Leibovici quite rightly wrote in comments, the exact solving involves the Airy functions. 
$$y(x)=c_1\text{Ai}\left(\sqrt[3]{\frac{a}{\epsilon^2}}\:x\right)+c_2\text{Bi}\left(\sqrt[3]{\frac{a}{\epsilon^2}}\:x\right)$$
Ai and Bi are the Airy functions http://mathworld.wolfram.com/AiryFunctions.html
The condition $y(\infty)=0$ is satisfied any $c_1$ and $c_2$, which shows that the condition $y(\infty)=0$ is not sufficient to fully determine a solution.
On physicists viewpoint, defining the properties of a boundary layer is certainly a way to introduce some boundary conditions in the range around $x=0$. Supposing that $y(x)\simeq y_0+y'_0 x$ in a layer of a given thickness $\delta$ allows to compute $c_1$ and $c_2$ as a function of $y_0$ , $y'_0$ and $\delta$, thanks to the series expansion of the Airy functions. Then the asymptotic expansion of the Airy functions will leads to the coefficient $C_2$ in the asymptotic approximation $(1)$. 
http://functions.wolfram.com/Bessel-TypeFunctions/AiryAi/06/
http://functions.wolfram.com/Bessel-TypeFunctions/AiryBi/06/
A: I found the example from my text book. I will add the details here: 


*

*First we see that in WKB approximation above
$$y(x) \sim \exp\left\{\frac{1}{\epsilon}S_0 + S_1\right\} = C_2 (ax)^{1/4} \exp\left[-\frac{1}{\epsilon}\sqrt{a} \frac{2}{3} x^{3/2} \right], \quad \epsilon \rightarrow 0^+,$$
we need to have 
$\epsilon S_2 << 1$ because generally we have $y(x) \sim \exp\left\{\frac{1}{\epsilon}S_0 + S_1 + \epsilon S_2 + \cdots \right\}$.
Solving for $S_2$ explicitly, we get $S_2 \sim C x^{-3/2}$, thus 
$\epsilon S_2 <<1 \Rightarrow \epsilon^{2/3} << x,$ and we have
$$\color{blue}{y_{one}(x) \sim C_2 (ax)^{1/4} \exp\left[-\frac{1}{\epsilon}\sqrt{a} \frac{2}{3} x^{3/2} \right] \quad \epsilon \rightarrow 0^+ \text{ on } x>> \epsilon^{2/3}}.$$

*Next, when $x$ is near zero, we let $t = (a/\epsilon^2)^{1/3}x$, so our ODE $\epsilon^2 y''(x) = axy(x)$ can be rewritten as $y''(t) = ty(t)$, and this is the Airy equation. And we know the solution is 
$$y(t) = D_1 \text{Ai}(t) + D_2 \text{Bi}(t)$$ 
thus 
$$\color{orange}{y_{two}(x) = D_1 \text{Ai}((a/\epsilon^2)^{1/3}x) + D_2 \text{Bi}((a/\epsilon^2)^{1/3}x)}.$$

*Finally would like to match them on $\epsilon^{2/3} <<x<<1$. First we see $y_{one}$ and $y_{two}$ are very different, so we will approximate $y_{two}$. From a irregular singular point  (ISP) analysis of the Airy equation at infinity, we will get 
$$\text{Ai}(t) \sim \frac{1}{2\sqrt{\pi}} t^{1/4} \exp\left[-\frac{2}{3}t^{3/2}\right], \quad t \rightarrow \infty $$
$$\text{Bi}(t) \sim \frac{1}{\sqrt{\pi}} t^{1/4} \exp\left[\frac{2}{3}t^{3/2}\right], \quad t \rightarrow \infty,$$
(the two constants $\frac{1}{2\sqrt{\pi}}$ and $\frac{1}{\sqrt{\pi}}$ are a bit more delicate and can not be obtained from the ISP analysis, but we will assume we have this fact.) So
$$\color{orange}{y_{two}(x) \sim \frac{D_1}{2\sqrt{\pi}} ((a/\epsilon^2)^{1/3}x)^{1/4} \exp\left[-\frac{2}{3}((a/\epsilon^2)^{1/3}x)^{3/2}\right] + \\
\frac{D_2}{\sqrt{\pi}} ((a/\epsilon^2)^{1/3}x)^{1/4} \exp\left[\frac{2}{3}((a/\epsilon^2)^{1/3}x)^{3/2}\right] \text{ when } x>>\epsilon^{2/3} \text{ and } \epsilon \rightarrow 0^+}. $$


After matching with $y_{one}$, we get $D_2 = 0$, $D_1 = 2\sqrt{\pi}(a\epsilon)^{-1/6} C_2$.
Now getting a uniform formula is non trivial, it follows from Langer's formula we have 
$$y_{uniform} (x) \sim C_2 2\sqrt{\pi}(a\epsilon)^{-1/6} \text{Ai}((a/\epsilon^2)^{1/3}x).$$
