I have to work out the integral of
$$ I(n):=\int_0^{\infty} u^n e^{-u} du $$
Somehow, the answer goes to
$$ I(n) = nI(n - 1)$$
and then using the Gamma function, this gives $I(n) = n!$
What I do is this:
$$ I(n) = \int_0^{\infty} u^n e^{-u} du $$
Integrating by parts gives
$$ I(n) = -u^ne^{-u} + n \int u^{n - 1}e^{-u} $$
Clearly the stuff in the last bit of the integral is now $I(n - 1)$, but I don't see how using the limits gives you the answer. I get this
$$ I(n) = \left( \frac{-(\infty)^n}{e^{\infty}} + nI(n - 1) \right) - \left( \frac{-(0)^n}{e^{0}} + nI(n - 1) \right) $$
As exponential is "better" than powers, or whatever its called, I get
$$ I(n) = (0 + I(n - 1)) + ( 0 + nI(n - 1)) = 2nI(n - 1)$$
Does the constant just not matter in this case?
Also, I do I use the Gamma function from here? How do I notice that it comes into play? Nothing tells me that $\Gamma(n) = (n - 1)!$, or does it?
You're doing something very problematic which is "evaluating at $\infty$". Improper integrals are to be evaluated as limits. So, you're looking at
$$I(n)=\lim_{m\to\infty}\int_0^m x^n e^{-x} dx=\int_0^\infty x^n e^{-x}dx$$
You're integrating by parts by writting an extra $nI(n-1)$ term (careful), so you ought to be writting
$$I(n) = \left(\lim_{m\to\infty} \frac{-m^n}{e^{m}} - \frac{-0^n}{e^{0}}\right) + nI(n - 1) $$
$$\underbrace {I(n)}_{udv} = \underbrace {\left( {\mathop {\lim }\limits_{m \to \infty } {{ - {m^n}} \over {{e^m}}} - {{ - {0^n}} \over {{e^0}}}} \right)}_{uv} + \underbrace {nI(n - 1)}_{vdu}$$
So that
$$I(n) = -\lim_{m\to\infty} \frac{m^n}{e^{m}} + nI(n - 1) $$
Now, can you evaluate
$$\lim_{m\to\infty} \frac{m^n}{e^{m}} ?$$
To make simpler, we could even do as follows. Set $$G(m,n)=\int_0^m x^n e^{-x}dx$$
Now, we integrate by parts
$$\displaylines{ G(m,n) = \int_0^m {{x^n}} {e^{ - x}}dx \cr = \left. { - {x^n}{e^{ - x}}} \right|_0^m + n\int_0^m {{x^{n - 1}}} {e^{ - x}}dx \cr = - {m^n}{e^{ - m}} + {0^n}{e^{ - 0}} + n\int_0^m {{x^{n - 1}}} {e^{ - x}}dx \cr = - {m^n}{e^{ - m}} + nG(m,n - 1) \cr} $$
Thus, taking $m\to \infty$,
$$\displaylines{ \Gamma \left( n+1 \right) = \mathop {\lim }\limits_{m \to \infty } G(m,n) \cr = \int_0^\infty {{x^n}} {e^{ - x}}dx \cr = - \mathop {\lim }\limits_{m \to \infty } {m^n}{e^{ - m}} + n\mathop {\lim }\limits_{m \to \infty } G(m,n - 1) \cr = 0 + n\Gamma \left( {n } \right) \cr = n\Gamma \left( {n} \right) \cr} $$
Without even considering the $\Gamma$ function, we get the desired recursion that implies $I(n)=n!$, or $\Gamma(n+1)=n!$. Note we just use $\Gamma$ because we know it "exists" but it'd be senseless to use the fact that $\Gamma(n+1)=n!$ to "solve" this, since we're precisely trying to figure out what $\Gamma(n+1)$ evaluates to.