$\renewcommand{\Re}{\operatorname{Re}}$ $\renewcommand{\Im}{\operatorname{Im}}$ $\newcommand{\bR}{\mathbb{R}}$ $\newcommand{\bC}{\mathbb{C}}$ $\newcommand{\bZ}{\mathbb{Z}}$ $\newcommand{\const}{\operatorname{const}}$
Consider now the same problem as before but with perturbation $M$ being of higher order than $L$. Then even for $\varepsilon \ll 1$ perturbation is not really small. In particular, for perturbed operator we need to have more boundary conditions than for unperturbed one and the different scenarios are possible. We consider only one scenario here when $L$ and $M$ are both positive: $Lu=\alpha(x)u$, $Mu=-(\beta (x)u')'$ with positive smooth $\alpha, \beta$: \begin{align} &L_\varepsilon u:=\alpha u -\varepsilon^2 (\beta (x)u')'=f\qquad 0\le x\le a, \label{eq-4.2.1}\\ &u(0)=g_1, \label{eq-4.2.2}\\ &u(a)=g_2. \label{eq-4.2.3} \end{align} Then boundary conditions (\ref{eq-4.2.2}) and (\ref{eq-4.2.3}) are important for $L_\varepsilon$ with $\varepsilon\ne 0$ but should be ignored for $L_0$. Recall that we assume \begin{equation} \alpha(x) >0, \qquad \beta (x) >0. \label{eq-4.2.4} \end{equation} This assumption guarantees that problem (\ref{eq-4.2.1})--(\ref{eq-4.2.3}) is well-posed.
Observe that we can construct by methods of the previous Section 4.1 \begin{equation} U_\varepsilon \sim \sum_{n\ge 0} U_n \varepsilon^n \label{eq-4.2.5} \end{equation} satisfying (\ref{eq-4.2.1}): \begin{align} &U_0=\alpha^{-1}f, \label{eq-4.2.6}\\ &U_{2m+2} = (\beta (\alpha^{-1}U_{2m})')' \label{eq-4.2.7}\\ &U_{2m+1}=0. \label{eq-4.2.8} \end{align} Then plugging $u= U_\varepsilon +v$ we get the same problem with $f=O(\varepsilon^\infty)$. Then actually $g_i$ will depend on $\varepsilon$ but it does not matter.
Remark 1. If $f=O(\varepsilon)$ then using Real Analysis one can prove easily that solution of the problem $L_\varepsilon w=0$, $w(0)=w(a)=0$ is $O(\varepsilon^\infty)$. Then plugging $U_\varepsilon:=U_\varepsilon+w$ we arrive to the same problem with $f=0$ exactly.
Remark 2. If $f=0$ then $u$ can reach positive maximum or negative minimum only on the ends of the interval $[a,b]$. Indeed, if $u$ reaches positive maximum in $c:0< c < a$, then $u(c) > 0$, $u'(c)=0$, $u''(c) \le 0$ and $L_\varepsilon u>0$.
So, $\min (g_1,g_2)\le u_\varepsilon \le \max (g_1,g_2)$.
Let us start from toy-model.
Example 1. Let $\alpha,\beta$ be constant. Then solution to (\ref{eq-4.2.1}) with $f=0$ is $u= C_1 e^{-\sigma x/\varepsilon} + C_2 e^{-\sigma (x-a)/\varepsilon}$ with $\sigma= (\alpha/\beta)^{-\frac{1}{2}}$ and plugging into (\ref{eq-4.2.2}) and (\ref{eq-4.2.3}) we can find $C_i=g_i + O(e^{-\sigma a/\varepsilon})$ and \begin{equation*} u= g_1 e^{-\sigma x/\varepsilon} + g_2e^{-\sigma (x-a)/\varepsilon}. \end{equation*}
We see that the ends should be treated separately and $u_\varepsilon $ is negligible outside of the boundary layers $\{x \gg \varepsilon \}$ and $\{a-x\gg \varepsilon\}$.
So we concentrate on the left end $x=0$ and we look at the solution in the form \begin{equation} u_\varepsilon = w e^{-\phi(x)\varepsilon^{-1}}. \label{eq-4.2.9} \end{equation} Then plugging into (\ref{eq-4.2.1}) and considering first only terms with $\varepsilon^0$ we get equation $\beta(\phi')^2 =\alpha$ and then \begin{equation} \phi(x) = \int_0^x (\alpha/\beta)^{\frac{1}{2}}\,dx. \label{eq-4.2.10} \end{equation} We take here $\phi(0)=0$ and $\phi(x)\sim \sigma x$ as $x>0$ with \begin{equation} \sigma=(\alpha(0)/\beta(0))^{\frac{1}{2}}. \label{eq-4.2.11} \end{equation} After $\phi(x)$ is fixed we consider all terms we get \begin{equation} \beta\varepsilon\Bigl( Pw+ \varepsilon Qw)\Bigr)e^{-\phi\varepsilon^{-1}}=0 \label{eq-4.2.12} \end{equation} with \begin{align} &Pw:=2\phi'w'+ (\phi''-\phi'\beta'\beta^{-1})w, \label{eq-4.2.13}\\ &Qw:=-w''+\beta'\beta^{-1}w'. \label{eq-4.2.14} \end{align} We can find asymptotic solution to this equation satisfying boundary condition \begin{equation} w|_{x=0}= g_1\sim \sum_{n\ge 0} g_{1n}\varepsilon ^n \label{eq-4.2.15} \end{equation} in the form \begin{equation} w\sim \sum_{n\ge 0} w_{n}\varepsilon ^n \label{eq-4.2.16} \end{equation} solving \begin{align} &Pw_n + Qw_{n-1}=0, \label{eq-4.2.17}\\ &w_n(0)=g_{1n}. \label{eq-4.2.18} \end{align}
One can consider different boundary condition, f.e. \begin{equation} \bigl(A \varepsilon u'+(B+C\varepsilon) u\bigr)|_{x=0}=0 \label{eq-4.2.19} \end{equation} Then (\ref{eq-4.2.18}) should be replaced by \begin{equation} (-A\sigma +B)w_n(0)= g_{1n}- Aw'_{n-1}(0)-Cw_{n-1}(0). \label{eq-4.2.20} \end{equation} Here $(-A\sigma +B)\ne 0$ is required but if it is not the case the trouble is much deeper than in our construction.