Consider the ODE $$\epsilon^2 y'' + \epsilon x y' - y = -1, \; y(0) = 0, \; y(1) = 3$$
I've seen the WKB method applied to homogeneous (linear) ODEs, but here we have the $-1$ term. I could perhaps do a change of variable writing $z = -y + 1$, but what if instead I consider the homogeneous case: $\epsilon^2 y'' + \epsilon x y' - y = 0$? Once I get the corresponding WKB solution how then do we obtain the WKB particular solution?
Regarding turning points, does $x = 0$ count as a turning point?
Nothing exciting happens. The WKB method provides approximations to basis solutions of the homogeneous ODE. You need to add a particular solution, here $y_p=1$, to the linear combinations of the basis solution to get the inhomogeneous solution.
Just applying the method, without special consideration for the second term at $x=0$, gives for the approximation of an exponential basis solution $y=\exp(S/ϵ)$, $S'=s=s_0+ϵs_1+...$, \begin{align} 0&=s^2+xs+ϵs'-1\\ (2s_0+x)^2&=4+x^2\\ 2s_0s_1+xs_1&=-s_0'=-\frac{x}{2s_0+x} \end{align} so in the first two terms $s_0=\frac12(-x\pm\sqrt{4+x^2})$ and $s_1=-\frac{x}{4+x^2}$, which integrates to give \begin{align} S_0&=-\frac14x^2\pm\frac14x\sqrt{4+x^2}+\ln(\sqrt{4+x^2}\pm x) \\ S_1&=-\frac12\ln(4+x^2) \end{align} The variant of $S_0$ that gives a bounded solution for $x>0$ is the one with the negative sign. However, as there is no outer solution from these basis solutions, there will also be a boundary layer at $x=1$. The full approximation is thus $$ y(x)=1-\frac{y_-(x)}{y_-(0)}+2\frac{y_+(x)}{y_+(1)}, $$ where $$ y_\pm(x)=\frac1{\sqrt{4+x^2}}\exp\left(\frac{\mp\frac{x}4(\sqrt{4+x^2}\pm x)+\ln(\sqrt{4+x^2}\pm x)}{ϵ}\right) $$