Consider \begin{equation} \varepsilon \frac{dy}{dx} = Q(x)y + R(x) \end{equation} where $\varepsilon$ is a small parameter. Can one apply a WKB method to find an asymptotic expansion for the solution?
I expect to obtain \begin{equation} y(x) \sim - \sum_{n = 0}^\infty \left(\frac{\varepsilon}{Q(x)} \frac{d}{dx}\right)^n \frac{R(x)}{Q(x)} \end{equation} (I am aware that this is a divergent sum), but I have been unsuccessful in recovering this form by using WKB methods.
On
This answer is to address a question in the comments of the other answer on whether the exponentially diverging terms drop out in regions where $\operatorname{Re} Q(x) < 0$.
We'll study the particular case of $Q(x) = -\sin x$ and $R(x) = 1$ over the interval $2\pi < x < 3\pi$. Note that we have $Q(x) < 0$ for all $x$ in question.
The general solution in $(1)$ becomes
$$ y(x) = y(0) e^{(\cos x - 1)/\epsilon} + \epsilon^{-1} e^{(\cos x - 1)/\epsilon} \int_0^x e^{(1-\cos \zeta)/\epsilon}\,d\zeta. \tag{2} $$
Let $\lambda = 1/\epsilon$, so that the integral in $(2)$ is
$$ I(x,\lambda) = \int_0^x e^{\lambda(1-\cos\zeta)}\,d\zeta. $$
We will obtain the asymptotic behavior of this integral as $\lambda \to \infty$ using the Laplace method.
When $2\pi < x < 3\pi$, the maximum value of $1-\cos\zeta$ on the interval $0 < \zeta < x$ is attained only when $\zeta = \pi$. Near there we have
$$ 1 - \cos\zeta = 2 - \tfrac{1}{2} (\zeta - \pi)^2 + O\left((\zeta-\pi)^3\right), $$
so by the Laplace method we obtain the asymptotic
$$ I(x,\lambda) \sim \int_{-\infty}^{\infty} \exp\left[\lambda(2 - \tfrac{1}{2}\zeta^2)\right]\,d\zeta = e^{2\lambda} \sqrt{\frac{2\pi}{\lambda}} $$
as $\lambda \to \infty$. Replacing $\lambda$ by $1/\epsilon$ we thus obtain an asymptotic for the term in $(2)$ corresponding to the particular solution:
$$ \epsilon^{-1} e^{(\cos x - 1)/\epsilon} \int_0^x e^{(1-\cos \zeta)/\epsilon}\,d\zeta \sim e^{(\cos x + 1)/\epsilon} \sqrt{\frac{2\pi}{\epsilon}} $$
as $\epsilon \to 0^+$.
So, while the homogeneous solution $y(0) e^{(\cos x - 1)/\epsilon}$ decays exponentially (since $\cos x - 1 < 0$ for $2\pi < x < 3\pi$), the remaining part of the solution still contains terms which diverge exponentially as $\epsilon \to 0^+$.
In general the differential equation is solved by
$$ y(x) = \exp\left(\frac{1}{\epsilon} \int_0^x Q(\xi)\,d\xi\right)\left[y(0) + \frac{1}{\epsilon}\int_0^x \exp\left(-\frac{1}{\epsilon} \int_0^\zeta Q(\xi)\,d\xi\right) R(\zeta)\,d\zeta \right]. \tag{1} $$
Unless we require $y(0) = 0$ then we will always have a term of the form $y(0)\exp\left(\frac{1}{\epsilon} \int_0^x Q(\xi)\,d\xi\right)$, and if the integral here is positive for some $x$ then this will not behave like the asymptotic series you propose as $\epsilon \to 0^+$.
The integral in the other term is of Laplace-type, and the largest/smallest contributions of its integrand come from neighborhoods of the roots of $Q(\zeta) = 0$ (assuming there are any).
Let's investigate the example $Q(x) = 2x$, $R(x) = 1$. In this case $(1)$ becomes
$$ y(x) = e^{x^2/\epsilon} \left[y(0) + \frac{1}{2}\sqrt{\frac{\pi}{\epsilon}} \operatorname{erf}\left(\frac{x}{\sqrt{\epsilon}}\right)\right]. $$
Assuming $x > 0$ we have
$$ \operatorname{erf}\left(\frac{x}{\sqrt{\epsilon}}\right) = 1 - \sqrt{\frac{\epsilon}{\pi}} \frac{e^{-x^2/\epsilon}}{x} + \frac{\epsilon^{3/2}}{2\sqrt{\pi}} \frac{e^{-x^2/\epsilon}}{x^3} + O\left(\epsilon^{5/2} e^{-x^2/\epsilon}\right) $$
as $\epsilon \to 0^+$, so
$$ y(x) = \frac{1}{2}\sqrt{\frac{\pi}{\epsilon}} e^{x^2/\epsilon} + y(0) e^{x^2/\epsilon} - \frac{1}{2x} + \frac{\epsilon}{4x^3} + O(\epsilon^2). $$
Notice that your expansion doesn't pick up either of the exponentially diverging terms but does agree with the part $1/2x + \epsilon/4x^3 + \cdots$. If you know you can ignore any diverging terms for your application then it looks like your series is the way to go. Otherwise it doesn't give you enough information.
As an example of how one might obtain the asymptotic behavior of the exponential terms we'll again take $Q(x) = 2x$ but this time let $R(x)$ be an arbitrary function which is analytic at $x = 0$. Then $(1)$ becomes
$$ y(x) = e^{x^2/\epsilon}\left(y(0) + \frac{1}{\epsilon}\int_0^x e^{-\zeta^2/\epsilon} R(\zeta)\,d\zeta \right). \tag{2} $$
Assuming $x > 0$ is fixed we can appeal to Watson's lemma and obtain an asymptotic expansion for the integral here by expanding $R(\zeta)$ in Taylor series around $\zeta = 0$. Using the fact that
$$ \int_0^\infty e^{-\zeta^2/\epsilon} \zeta^k\,d\zeta = \frac{1}{2} \Gamma\left(\frac{k+1}{2}\right) \epsilon^{(k+1)/2} $$
we get
$$ y(x) \approx e^{x^2/\epsilon}\left(y(0) + \frac{1}{2}\sum_{k=0}^{\infty} \frac{R^{(k)}(0) \Gamma\left(\frac{k+1}{2}\right)}{k!} \epsilon^{(k-1)/2}\right) $$
as $\epsilon \to 0^+$.