According to Sturm-Liouville's theory, any second order linear ordinary differential operator $L=P(x)\dfrac{d^2}{dx^2}+Q(x)\dfrac{d}{dx}+R(x)$ can be reduced to the Sturm–Liouville form, $\mathcal{L}=\dfrac{d}{d x}\left[p(x) \dfrac{d}{d x}\right]+q(x)$, using an integrating factor $\mu(x)$ given by
$$ \mu(x)=\frac{1}{P(x)} \exp \left(\int \frac{Q(x)}{P(x)} d x\right) $$
Then, $p(x)=\mu(x)P(x)$, $q(x)=\mu(x)R(x)$. This allows us to find a set of eigenfunctions $y_n$ for the Sturm-Liouville problem, $\mathcal{L}y_n=\lambda_n y_n$, and compute the solution of a problem $\mathcal{L}y=f(x)$ spanning it in the eigenfunction basis.
However, if we wanted the solution of our original problem $Ly=F(x)$, with $F(x)=f(x)/\mu(x)$, how could the solutions of this problem be obtained from the associated Sturm-Liouville problem?
For example. If we have $xe^xy''+e^xy'+\dfrac{e^x}{x}y=\dfrac{e^x}{x}$, with $x\in[1,e]$ and $y(1)=y(e)=0$, it can be reduced to the S-L form using $\mu(x)=e^{-x}$. Then, the normalized eigenfunctions of the S-L problem would be
$$\phi_n(x)=\frac{1}{\sqrt{2}}\sin(n\pi\ln x), \quad n=1,2,\ldots$$
And, after some calculations, the solution of the inhomogeneous S-L problem would turn out to be
$$ y(x)=\sum_{n=1}^{\infty} \frac{2}{n \pi} \frac{\left[(-1)^{n}-1\right]}{n^{2} \pi^{2}-1} \sin (n \pi \ln (x)) $$
How could we recover the solutions of the original problem from these solutions of the S-L problem?