Solution (source code)

= Solution

Assume first that $n\geq m$. Substituting <Rodrigues' formula> for $L_n$ gives
$$
I_{nm}=\frac1{n!}\int_0^\infty L_m(x)\frac{d^n}{dx^n}(x^ne^{-x})\,dx
=\frac{(-1)^n}{n!}\int_0^\infty L_m^{(n)}(x)x^ne^{-x}\,dx.
$$
Here <integration by parts> was performed $n$ times. Each boundary term vanishes: at infinity a <polynomial> times $e^{-x}$ tends to zero; at zero every <derivative> of $x^ne^{-x}$ of order less than $n$ vanishes. If $n>m$, the <derivative> $L_m^{(n)}$ is zero, proving <orthogonality>. The case $m>n$ follows by exchanging the two factors.

For $m=n$, the leading coefficient found above gives $L_n^{(n)}=(-1)^n$. Therefore
$$
I_{nn}=\frac1{n!}\int_0^\infty x^ne^{-x}\,dx=\frac{\Gamma(n+1)}{n!}=1.
$$
The same conclusion includes $n=0$ directly. Thus the full <Laguerre orthonormality> relation is
$$
\boxed{\int_0^\infty L_n(x)L_m(x)e^{-x}\,dx=\delta_{nm}.}
$$