WKB for a boundary value problem with two layers The equation $\epsilon y''-x^4y'-y=0,\ y(0)=y(1)=1$, can be solved by boundary layer analysis and turns out it has two layers of size $O(\epsilon) $ and $O(\sqrt{\epsilon})$ at $1,0$ accordingly. Is possible to solve this with WKB method?
How do I deal with the fact that the coefficient of $y',\ q(x)=-x^4$ is $q(0)=0$?

Boundary layer analysis:
outter expansion :$\epsilon\rightarrow0 \implies -x^4y'-y=0\implies y_o(x)=Ce^{\frac{1}{3x^3}}$ 
Since $q(x)\leq0$, we expect a boundary layer near $1$ . Substituting $X=\frac{(1-x)}{\delta(\epsilon)}\implies$
$\frac{\epsilon}{\delta^2(\epsilon)}Y''(X)+(1-\epsilon X)^4\frac{1}{\delta(\epsilon)}Y'(X)-Y(X)=0.$ Scaling $\delta(\epsilon)=O(\epsilon)$ as $\epsilon\rightarrow0\implies$
$Y''(X)+Y'(X)=0\implies Y_1(X)=A_0+B_0e^{-X}$. $y(1)=1\implies Y(0)=1\implies A_0+B_0=1$.
From $y_o(x)$ we can see that as $x\rightarrow0,\ y_o(x)\rightarrow\infty$ unless $C=0$ thus $y_o(x)=0$. 
Matching to outter solution $\lim_{X\rightarrow\infty}Y_1(X)=0\implies A_0=0\implies B_0=1$ and $Y_1(X)=e^{-X}.$
Then there is a boundary layer at $0$ obviously.
 Substituting $Z=\frac{x}{\delta(\epsilon)}\implies$ $\frac{\epsilon}{\delta^2(\epsilon)}Y''(Z)-(\delta(\epsilon)Z)^4\frac{1}{\delta(\epsilon)}Y'(Z)-Y(Z)=0$.
The only consistent scale is $\delta(\epsilon)=O(\sqrt{\epsilon})$.
$Y(Z)''-Y(Z)=0\implies Y_0(Z)=D_0e^Z+E_0e^{-Z}$. $y(0)=1\implies D_0+E_0=1$.
Again, matching to the outer solution $\lim_{Z\rightarrow\infty}Y_0(Z)=0$ is only possible if $D_0=0$. Thus $Y_0(Z)=e^{-Z}$, and finally
$y_{uniform}(x)=e^{\frac{x-1}{\epsilon}}+e^{-\frac{x}{\sqrt{\epsilon}}}$.
 A: The WKB variation that I am aware of for this class of problem begins by using $y=fg$ to remove the first derivative terms of $g$ by suitable choice of $f$. So here, you want
$$\epsilon (f'' g + 2f' g' + fg'')-x^4(f'g + g'f)-fg=0$$
and you want the $g'$ terms to go away. So you choose $2\epsilon f' - x^4 f = 0$ and then go for $f=\exp(x^5/(10\epsilon))$. Then the equation becomes
$$\epsilon g'' + \left ( 2x^3-1-\frac{x^8}{4\epsilon} \right ) g = 0.$$
Now with $g=\exp(\theta/\epsilon) h$ you have
$$4\epsilon^2 h''+(((8x^3-4+4\theta'')h+8\theta'h')\epsilon+(4\theta'^2-x^8) h=0.$$
Now the leading order equation is apparently $\theta'^2=x^8/4$. But if this is true all the time then for small enough $x$, $\theta'$ becomes much smaller than $\epsilon$, and so it is not appropriate to discard, for example, the $-4\epsilon$ term. In view of that, you might hope that it is consistent to keep that $-4\epsilon$ term only, i.e. to have
$$\theta'^2=x^8/4+\epsilon$$
so $\theta(x)=\int_0^x \sqrt{\epsilon+y^8/4} dy$ is a valid choice, resulting in the overall ansatz $y(x)=e^{\frac{\frac{x^5}{10}+\int_0^x \sqrt{\epsilon+y^8/4} dy}{\epsilon}} h(x)$. (Which solution you select for $\theta$ does not matter because it only affects the overall solution by a constant factor that can be absorbed into $h$.) So the argument of the exponential smoothly transitions from behaving like $x/\sqrt{\epsilon}$ to behaving like $x^5/(5\epsilon)$ as you pass through a layer of total thickness $O(\epsilon^{1/8})$. This correctly produces the boundary layer at $x=0$ assuming $h$ is regular there.
Accordingly you get
$$4\epsilon^2 h''+ 8\epsilon \sqrt{\epsilon+\frac{x^8}{4}} h' + \epsilon \left ( 8x^3+\frac{4x^7}{\sqrt{\epsilon+x^8/4}} \right ) h=0 \\
h(0)=1 \\
h(1)=e^{-\frac{1}{\epsilon} \left ( \frac{1}{10} + \int_0^1 \sqrt{\epsilon+y^8/4} dy \right )}$$
which is arguably nastier than what you started with, even though it does only have one boundary layer now. I guess the next step is to remove the $h'$ term by applying the same trick and then repeat. The result seems a lot nastier than the boundary layer analysis, though it retains a feature (the behavior in the $\epsilon^{1/2} \ll x \ll \epsilon^{1/8}$ region) which is lost in the boundary layer analysis.
When I tried direct WKB analysis on the original equation, it seemed to fail spectacularly, because I could not isolate terms only involving $y_0$ rather than any of its derivatives to determine $\theta$. Indeed if $y=\exp(\theta/\epsilon) h$ then $y'=\exp(\theta/\epsilon)((\theta'/\epsilon) h+h')$ and $y''=\exp(\theta/\epsilon)((\theta'/\epsilon)((\theta'/\epsilon) h+h')+(\theta''/\epsilon) h + (\theta'/\epsilon) h' + h'')$ so you get
$$\epsilon ((\theta'/\epsilon)((\theta'/\epsilon) h+h')+(\theta''/\epsilon) h + (\theta'/\epsilon) h' + h'') -x^4((\theta'/\epsilon) h+h')-h=0.$$
For non-small $x$ you can group up the terms of order $1/\epsilon$ here and get $\theta'^2=x^4\theta'$ which seems to lead to $\theta=\frac{x^5}{5}$. But for small $x$, the $x^4$ term is no longer dominant, and so another term has to be kept to determine $\theta$. But how do you pick which one? I couldn't figure it out, but maybe there is a way.
A: Let's try to hack through with the most simplistic approach, search for a basis function $y=e^S$, with the guiding idea that $|s'|\ll s^2$, $s=S'$. Then the equation becomes
$$
0=ϵ(s'+s^2)−x^4s−1.
$$
By the guiding idea, solve the quadratic equation $0=ϵs^2−x^4s−1$,
$$
s_1=\frac{x^4+\sqrt{x^8+4ϵ}}{2ϵ}\approx\frac{x^4}{ϵ}
$$
or
$$
s_2=\frac{-2}{x^4+\sqrt{x^8+4ϵ}}
=-\frac1{\sqrtϵ}+\frac{x^4+(\sqrt{x^8+4ϵ}-\sqrt{4ϵ})}{\sqrtϵ(x^4+\sqrt{x^8+4ϵ})}
=-\frac1{\sqrtϵ}+\frac{x^4}{\sqrtϵ}·\frac{1+\frac{x^4}{\sqrt{x^8+4ϵ}+\sqrt{4ϵ}}}{x^4+\sqrt{x^8+4ϵ}}
$$
Thus $y_1=\exp(S_1)\approx\exp(\frac{x^5-1}{5ϵ})$ has a boundary layer of width $1$ at $x=1$, while $S_2(x)=-\frac{x}{\sqrtϵ}-...$ gives a layer at $x=0$ of width $1/2$ with a second term that becomes important further inwards at about $x=\sqrt[8]ϵ$. For larger $x$ the full expression does not become zero, one gets $S_2(x)\approx c+x^{-3}/3$.
I do not see what a suitable value of $c$ or a matching formula could be, also the exact integral formula is not immediately helpful,
$$
S_2(x) = \frac{x}{ϵ}\,{}_2F_1\left(-\frac12, \frac18, \frac98, -\frac{x^8}{4ϵ}\right) - \frac{x^5}{10ϵ} + c
$$
