Elementary proof for $\lim\limits_{n \to\infty}\frac{n!e^n}{n^n} = +\infty$

Considering the ratio of consecutive elements of $$ a_n=\frac{n!e^n}{n^n}\tag{1} $$ we get $$ \begin{align} \frac{a_n}{a_{n-1}} &=\frac{(n-1)^{n-1}e}{n^{n-1}}\\ &=e\left(1-\frac1n\right)^{n-1}\tag{2} \end{align} $$ Taking the log of this ratio gives $$ \begin{align} \log\left(\frac{a_n}{a_{n-1}}\right) &=1+(n-1)\log\left(1-\frac1n\right)\\ &=1-(n-1)\left(\frac1n+\frac1{2n^2}+O\left(\frac1{n^3}\right)\right)\\ &=\frac1{2n}+O\left(\frac1{n^2}\right)\tag{3} \end{align} $$ Summing $(3)$ yields $$ \log(a_n)=\frac12H_n+C+O\left(\frac1n\right)\tag{4} $$ Since the harmonic series diverges, $(4)$ shows that $\log(a_n)$, and hence $a_n$, grows without bound.

One can use integrals. The work is a little long, but quite rewarding. Consider $$I_n=\int_0^{\pi /2}\sin^{n}tdt$$

Then $$I_{2n+1}\leqslant I_{2n}\leqslant I_{2n-1}$$

On the other hand integrating by parts gives

$${I_{2n + 1}} = \frac{{2n}}{{2n + 1}}{I_{2n - 1}}$$

It follows that $$1\leqslant \frac{{{I_{2n}}}}{{{I_{2n + 1}}}}\leqslant \frac{{{I_{2n - 1}}}}{{{I_{2n + 1}}}}$$ so $$\frac{{{I_{2n}}}}{{{I_{2n + 1}}}}\to 1$$

We proceed to evaluate that limit in a different way. Integrating by parts like before gives us $${I_{2n}} = \frac{{\left( {2n - 1} \right)!!}}{{\left( {2n} \right)!!}}\frac{\pi }{2}$$ $${I_{2n + 1}} = \frac{{\left( {2n } \right)!!}}{{\left( {2n+ 1} \right)!!}}$$

Thus, we find $$\frac{{{I_{2n + 1}}}}{{{I_{2n}}}} = \frac{{\left( {2n} \right)!{!^2}}}{{\left( {2n - 1} \right)!!\left( {2n + 1} \right)!!}}\frac{2}{\pi } \to 1$$ or $$\prod\limits_{k = 1}^\infty {\frac{{4{k^2}}}{{4{k^2} - 1}}} = \frac{\pi }{2}$$

This is the celebrated Wallis product for $\pi$. Now onto Stirling. Consider the limit $$A = \mathop {\lim }\limits_{n \to \infty } \frac{{n!{e^n}}}{{{n^{n + 1/2}}}}$$ Then one must have by $n\mapsto 2n$ that $$1 = \frac{A}{A} = \mathop {\lim }\limits_{n \to \infty } \frac{{n!{e^n}}}{{{n^{n + 1/2}}}}\frac{{{{\left( {2n} \right)}^{2n + 1/2}}}}{{\left( {2n} \right)!{e^{2n}}}} = \sqrt 2 \mathop {\lim }\limits_{n \to \infty } \frac{{{n^n}}}{{{e^n}n!}}\frac{{\left( {2n} \right)!!}}{{\left( {2n - 1} \right)!!}}$$

We now use Wallis product formula. Squaring, this means that $$1 = 2{\left( {\mathop {\lim }\limits_{n \to \infty } \frac{{{n^{n + 1/2}}}}{{{e^n}n!}}} \right)^2}\left( {\mathop {\lim }\limits_{n \to \infty } \frac{{2n + 1}}{n}} \right)\left( {\mathop {\lim }\limits_{n \to \infty } \frac{{\left( {2n} \right)!!\left( {2n} \right)!!}}{{\left( {2n - 1} \right)!!\left( {2n + 1} \right)!!}}} \right)$$ And we thus get Stirling's $$\frac{1}{{\sqrt {2\pi } }} = { {\mathop {\lim }\limits_{n \to \infty } \frac{{{n^{n + 1/2}}}}{{{e^n}n!}}} }$$

ADD Denote by $a_n$ our sequence above. Then $${a_n} > {a_{n + 1}}$$ is equivalent to $${\left( {1 + \frac{1}{n}} \right)^{n + 1/2}} > e$$ under rearranging. It remains to prove this.

Similar to robjohn's explanation, one can approach this using the so-called "trapezoidal rule" for estimating integrals.

Thus, consider the integral $\int_1^n \log(x) \; dx$. The trapezoidal rule estimate gives, over each interval $[m, m+1]$, the estimate

$$\int_m^{m+1} \log(x) \; dx \;\; = \;\; \frac1{2}(\log(m) + \log(m+1)) + e_m$$

where the error term $e_m$ is bounded above by a constant times $m^{-2}$. Summing from $m=1$ to $n-1$, we arrive at

$$x\log(x) - x\; |_1^n = \int_1^n \log(x)\; dx \;\; = \;\; \frac1{2}\log(n)+ \sum_{1}^{n-1} \log(n) + E_n$$

where the error term $E_n = \sum_{m=1}^{n-1} e_m$ is some constant $E = \sum_{m=1}^\infty e_m$ minus an amount that is bounded above by some constant $C$ times $n^{-1}$. We may exponentiate to get $n^n e^{-n}e^{-1} = e^{E_n}(\frac{n!}{\sqrt{n}})$, which may be massaged into the asymptotic formula

$$\frac{e^n n!}{n^n} = (1 + O(n^{-1}))e^{1-E}\sqrt{n}$$

where $E$ could be determined if we want, but this gets the job done. See Terry Tao's post here.