Solution (source code)

= Solution

For two real <eigenfunctions>, multiply the equation for $y_n$ by $y_m$ and the equation for $y_m$ by $y_n$, then subtract. The $q$ terms cancel. The resulting <Lagrange identity for a Sturm-Liouville operator> is
$$
\left[p(y_my_n'-y_ny_m')\right]'+(\lambda_n-\lambda_m)r y_my_n=0.
$$
Integrating first over $[\epsilon,R]$ and then taking the endpoint limits gives
$$
(\lambda_n-\lambda_m)\int_0^\infty r y_my_n\,dx=-\left[p(y_my_n'-y_ny_m')\right]_0^\infty.
$$
The needed assumptions are that the <improper integral> exists and that the displayed boundary expression tends to zero at both endpoints. Bounded $y_m,y_n$ and their first <derivatives> are sufficient because $p\to0$, but are stronger than necessary. Merely knowing that $p\to0$ does not control arbitrarily growing <eigenfunctions>. With the vanishing boundary form and distinct <eigenvalues>, division by $\lambda_n-\lambda_m$ proves <orthogonality of Sturm-Liouville eigenfunctions>:
$$
\boxed{\int_0^\infty r(x)y_m(x)y_n(x)\,dx=0\quad(m\ne n).}
$$
For complex <eigenfunctions> in a real <self-adjoint differential operator> problem, use $\overline{y_m}$ instead, giving the corresponding <weighted inner product>.

Multiplying the particular equation by $e^{-x}$ puts it in <Sturm-Liouville form>:
$$
(xe^{-x}y')'+\lambda e^{-x}y=0.
$$
Consequently
$$
\boxed{p(x)=xe^{-x},\qquad q(x)=0,\qquad r(x)=e^{-x}.}
$$
Its polynomial <eigenfunctions> are the <Laguerre polynomials>. Products of these <polynomials>, with or without one <derivative>, grow only algebraically at infinity. The exponential in $p$ therefore makes the boundary form vanish there; its factor $x$ makes it vanish at zero. Their <orthogonality> is thus
$$
\int_0^\infty e^{-x}y_m(x)y_n(x)\,dx=0\quad(m\ne n).
$$

For $\lambda=1$, substitution of $Ax+B$ gives $A+B=0$. For $\lambda=2$, substitution of $Ax^2+Bx+C$ gives $4A+B=0$ and $B+2C=0$. Choosing $y_n(0)=1$, which fixes the otherwise arbitrary nonzero multiplicative constants, yields
$$
\boxed{y_1=1-x,\qquad y_2=1-2x+\frac{x^2}{2}.}
$$
Their product is $1-3x+\tfrac52x^2-\tfrac12x^3$. Using the <Gamma function> identity $\int_0^\infty x^ke^{-x}dx=k!$, also obtainable by repeated <integration by parts>, verifies the required <orthogonality> directly:
$$
\int_0^\infty e^{-x}y_1y_2\,dx=1-3(1!)+\frac52(2!)-\frac12(3!)=0.
$$