I am working on a boundary layer problem for a second order linear ODE. A simpler problem which I think still illustrates the issue I am having is
$$\varepsilon y''-y'+y=0,y(0)=0,y(1)=1$$
where $\varepsilon > 0$ is a small parameter. This problem (unlike my actual problem) in fact admits an exact solution, namely
$$\frac{e^{\lambda_2 x}-e^{\lambda_1 x}}{e^{\lambda_2}-e^{\lambda_1}}$$
where
$$\lambda_1,\lambda_2=\frac{1 \pm \sqrt{1-4\varepsilon}}{2 \varepsilon}.$$
For small $\varepsilon > 0$, $\lambda_1$ is a large positive number, namely $\frac{1}{\varepsilon}+O(1)$, while $\lambda_2$ is a positive number of order $1$, namely $1+O(\varepsilon)$.
Upon sorting out minus signs, this means that the solution above grows very fast near $1$ and is very small far away from $1$. Accordingly, I've developed an "inner" boundary layer solution near $x=1$. This goes through the usual way by changing variables to $z=\frac{x-1}{\varepsilon}$, assuming that $\frac{d^2 y}{dz^2},\frac{dy}{dz}$ and $y$ are all $O(1)$ near $x=1$, and balancing the leading order terms. From this I get a leading order solution of $e^{\frac{x-1}{\varepsilon}}$ for $|x-1|=O(\varepsilon)$.
Now my problem arises on the rest of the interval. For the leading order outer solution, we consider $y'=y$, and get $y=Ce^x$. My problem is now essentially in identifying $C$. I can take it to be $0$ to match the left boundary condition, but this doesn't seem right, because the composite solution from this procedure is $e^{\frac{x-1}{\varepsilon}}$ over the whole interval, which fails to satisfy the left boundary condition (at least exactly). I can subtract off $e^{-\frac{1}{\varepsilon}}$ to force the left boundary condition to hold (which has a negligible impact on the error in the DE itself), but then the right boundary condition doesn't hold.
Should I be thinking about this differently? Do I need to go to higher order in order to get a reasonable result? It seems to me that I will, because a "reasonable" first order linear ODE beginning with value zero will always just stay at zero. So to start at zero and not stay there I need to consider a truly second order problem. But this is more difficult in my actual problem, and it seems that the method should be more or less universal.
The WKB expansion is $$ y(x)=\exp\left(\frac{1}{\delta}\sum_{n=0}^{\infty}\delta^nS_n(x)\right), $$ where $\delta=\epsilon^\alpha$. The derivative of $y$ is $$ y'(x)=\left[\frac{1}{\delta}\sum_{n=0}^{\infty}\delta^nS_n'(x)\right]y(x), $$ and the second derivative is, $$ y''(x)=\left[\frac{1}{\delta^2}\left(\sum_{n=0}^{\infty}\delta^nS_n'(x)\right)^2+\frac{1}{\delta}\sum_{n=0}^{\infty}\delta^nS_n''(x)\right]y(x). $$
So, substituting into the ODE and factoring out the exponentials, you get $$\left[\frac{\epsilon}{\delta^2}\left(\sum_{n=0}^{\infty}\delta^nS_n'(x)\right)^2+\frac{\epsilon}{\delta}\sum_{n=0}^{\infty}\delta^nS_n''(x)\right]-\left[\frac{1}{\delta}\sum_{n=0}^{\infty}\delta^nS_n'(x)\right]+1=0. $$ You now have to find a distinguished limit to find $\delta$. The two largest terms in the expansion are $S_0'(x)^2\epsilon/\delta^2$ and $S_0'(x)/\delta$. To balance these, choose $\delta=\epsilon$. (This is the same as the boundary layer scaling you get from normal boundary layer theory).
Now, pick out the leading order eikonal equation, at $O(\epsilon^{-1})$: $$ \left(S_0'(x)\right)^2-S_0'(x)=0,$$ and the first order transport equation, at $O(1)$: $$ 2S_0'(x)S_1'(x)+S_0''(x)-S_1'(x)+1=0.$$
The first equation gives two solutions, $S_0(x)=x$, and $S_0(x)=0$. The first solution is an inner solution, and the second is an outer solution. You don't need to add constants of integration. The leading order solution is $$y(x)=A\exp\left(\frac{x}{\epsilon}\right)+B, $$ and boundary conditions give $$ y(x)=\frac{1-\exp\left(\frac{x}{\epsilon}\right)}{1-\exp\left(\frac{1}{\epsilon}\right)}. $$
The solution to the transport equation is, for the $S_0(x)=0$ solution, $$ S_1(x)=x,$$ and for the $S_0(x)=x$ solution, $$ S_1(x)=-x.$$ This gives the two-term solution, $$ y(x)=C\exp\left(x\right)+D\exp\left(\frac{x}{\epsilon}-x\right). $$ Boundary conditions give the final two-term solution as $$ y(x)=\frac{\exp\left(x\right)-\exp\left(\frac{x}{\epsilon}-x\right)}{\exp\left(1\right)-\exp\left(\frac{1}{\epsilon}-1\right)}. $$ A nice thing about WKB theory is that it gives you the inner and outer solutions together, you don't need to do asymptotic matching.
Bender and Orszag's book has a great section about WKB theory.
Here's a picture of the one and two term solutions, along with a numerical solution, for $\epsilon=0.2$.