Solution (source code)

= Solution

Choose a complete orthonormal system of <eigenfunctions> $\phi_n$ for the <self-adjoint> problem. Expand $f=\sum_n f_n\phi_n$ and $y=\sum_n c_n\phi_n$, with
$$
f_n=\langle\phi_n,f\rangle=\int_a^b\phi_n^*(\xi)f(\xi)w(\xi)\,d\xi.
$$
Substitution into $Ly=f$ gives $\lambda_nc_n=f_n$. Since no <eigenvalue> is zero, the <Sturm-Liouville eigenfunction expansion> is
$$
\boxed{y(x)=\sum_n\frac{\phi_n(x)}{\lambda_n}\int_a^b\phi_n^*(\xi)f(\xi)w(\xi)\,d\xi.}
$$
For <eigenfunctions> not normalized to one, divide each term additionally by $\langle\phi_n,\phi_n\rangle$. The expansion is understood in the weighted function space and, with sufficient data regularity, in the operator domain. A zero <eigenvalue> would instead require its corresponding source coefficient to vanish and would leave an arbitrary homogeneous component.