Solution (source code)

= Solution

For $n>1$, <integration by parts> with $u=x^{n-1}$ and $dv=e^{-x}\,dx$ yields
$$
I_n=[-x^{n-1}e^{-x}]_0^\infty+(n-1)\int_0^\infty x^{n-2}e^{-x}\,dx=(n-1)I_{n-1}.
$$
The boundary term vanishes at zero because $n-1>0$, and at infinity because exponential decay dominates every fixed power. The initial <improper integral> is $I_1=\int_0^\infty e^{-x}\,dx=1$. Iterating the <recurrence relation> therefore gives
$$
\boxed{I_n=(n-1)(n-2)\cdots1\,I_1=(n-1)!,\qquad n\ge1.}
$$
This is the positive-integer case of the <Gamma function> identity $\Gamma(n)=(n-1)!$.