Construct leading order approximation without given specific function The Question:
$$\varepsilon y''+f(x)y'+y=0 \qquad y(-1)=0 \qquad y(1)=1$$
where $0<\varepsilon \ll 1$ and $f$ is a given smooth function that is strictly positive with $f(1)=f(-1)=1$.
(i) Determine the location of the boundary layer
(ii) Obtain leading order outer and inner solutions

My Attempt:
(i) By first assuming that $y(x,\varepsilon) \sim y_0(x) + \varepsilon y_1(x)+\cdots$ is of order $O(1)$, I plugged this into the equation to obtain
$$f(x)y_0'(x)+y_0(x)=0 \qquad y_0(-1)=0 \qquad y_0(1)=1$$
at the leading order. Solving this gives
$$y_0(x) = A\exp\biggl(-\int \frac{dx}{f(x)}\biggl)$$
for some constant $A$. If we try to satisfy $y_0(-1)=0$ here, we get $A=0$ which is a contradiction, for then $y$ would not be of order $O(1)$.
It follows that the boundary layer is at $x=-1$.
(ii) So I have already found $y_0$ above, which is indeed the leading order term to the outer solution.
But how do I find the constant $A$ if I don't know what $f$ is?
How do I use the fact that $f(-1)=f(1)=1$?
 A: As an alternative to contemplating inner and outer solutions, the WKB approximation method directs you to consider $y=\exp(S/\delta)$ with accordingly
$$
ε(δS''+S'^2)+δfS'+δ^2=0
$$
By the assumption $f>0$, and thus $\min_{x\in[0,1]} f(x)\gg \max(ε,δ)$, the leading terms in magnitude are $εS'^2+δfS'$ and they are in magnitude balance for $ε=δ$. Then the perturbation problem to solve is
$$
S'^2+fS'=-ε(S''+1).
$$ 
The series expansion $S=S_0+εS_1+ε^2S_2$ results in (letting $s_k=S_k'$) two solution for the zeroth order equation $s_0^2+fs_0=0$ that lead to two expansions giving approximations to two independent basis solutions.
\begin{align}
s_0^2+fs_0&=0&:&&s_0&=0&\text{or}&& s_0&=-f(x)&
\\
(2s_0+f)s_1&=-(s_0'+1)&:&&s_1&=-\frac1f&|&&
s_1&=\frac{-f'+1}f
\\
s_1^2+(2s_0+f)s_2&=-s_1'&:&&s_2&=-\frac{1+f'}{f^3}&|&&s_2&=\frac{-f''f+(-2f'+1)(-f'+1)}{f^3}
\end{align}
etc.
Let $F$ be an anti-derivative of $f$, $G$ of $1/f$, both are monotonously increasing. Select the integration constants so that $F(-1)=0=G(-1)$. Then the two basis solutions using $s_0$ and $s_1$ give the approximation
$$
y(x)\approx Ae^{-G(x)}+\frac{B}{f(x)}e^{-\frac1εF(x)+G(x)}
$$
The boundary conditions imply
\begin{align}
0&=A+B&\implies A&=-B\\
1&=-Be^{-G(1)}+Be^{-\frac1εF(1)+G(1)}&\implies B&=-e^{G(1)}+\text{very small terms},
\end{align}
the last because $\exp(-\frac1εF(1))$ is $O( ε^k)$ for any order $k$. In total
$$
y(x)\approx e^{G(1)-G(x)}-\frac1{f(x)}e^{-\frac1εF(x)+G(1)+G(x)}.
$$
The first part corresponds to the outer solution, the second the inner solution for the boundary layer at $x=-1$. It is a little more complex than the inner solution $\exp(-\frac1εx+G(1))$ one would get with the direct method.
