Here is a proof based on the discussion under David E Speyer's answer. It uses the same method also to sharpen the proof of the lower bound.
We will give a simple proof for the bound $$ (n/e)^n\sqrt{2\pi n} \le n! \le (n/e)^n(\sqrt{2\pi n}+1). $$
We use the integral representation of $n!$: \begin{align*} n! &= \int_0^\infty x^n e^{-x} %= \int_0^\infty e^{n \log(x) - x} = n\int_0^\infty e^{n \log(nx) - nx} dx %= n n^{n}\int_0^\infty e^{n (\log(x) - x)} dx = n (n/e)^{n}\int_{-1}^\infty e^{n (\log(1+x) - x)} dx. \end{align*} Note that $\log(1+x)-x$ is maximized at $x=0$. Inspired by the series expansion $\log(1+x)-x=-x^2/2 + O(x^3)$ we make a change of variables, $\log(1+x)-x=-z^2/2$, or: $$ z = \text{sign}(x)\sqrt{2(x-\log(1+x))}. $$ The differentials turn out to be pretty simple: $$ d(-z^2/2)= -zdz = d(\log(1+x)-x) = \frac{dx}{1+x} - dx = -\frac{x}{1+x}dx, $$ which after adjusting the integration limits yield: $$ n! = (n/e)^{n} n\int_{-\infty}^\infty e^{-nz^2/2} (z+z/x) dz. $$ We now need the bounds \begin{align} 1+2z/3 &\le z+z/x \le \max\{1, 1+z\} \quad\text{for all }x\in\mathbb R. &&(0) \end{align} Since $\int_{-\infty}^\infty e^{-nz^2/2}dz=\sqrt{2\pi/n}$ and $\int_{0}^\infty e^{-nz^2/2} z\, dz=1/n$, this immediately gives us $$ \sqrt{2\pi n} \le n\int_{-\infty}^\infty e^{-nz^2/2} (z+z/x) dz \le \sqrt{2\pi n} + 1. $$ Which is what we want.
It remains to show eq. (0), which we will do using the following standard logarithmic inequalities: \begin{align} \log(1+x) &\ge x(2+x)/(2+2x) \quad &\text{when } x<0 &&(1) \\ \log(1+x) &\ge x-x^2/2 \quad &\text{when } x>0 &&(2) \\ \log(1+x) &\le \frac{x(6+x)}{6+4x} \quad &\text{when } x>-1 &&(3) \end{align} We note $z(1+1/x)$ is non-negative, since $1+1/x<0$ implies $\mathrm{sign}(x)=-1$. Thus it suffices to show $z^2(1+1/x)^2\le 1$, which follows using eq. (1): $$ z^2 = 2(x-\log(1+x)) \le \frac{x^2}{1+x} \le \frac{x^2}{(1+x)^2} = (1+1/x)^{-2}, $$ For $x>0$ we can use eq. (2) to bound $ z = \sqrt{2(x-\log(1+x))} \le \sqrt{2(x^2/2)} = x, $ so $z/x\le 1$. This proves the upper bound of (0). For the lower bound, it suffices to show $z^2(1/3+1/x)^2 \ge 1.$ Using eq. (2) we get \begin{align*} z^2(1/3+1/x)^2 %&= 2(x - \log(1+x))(1/3+1/x)^2 &\ge 2(x - x(6+x)/(6+4x))(1/3+1/x)^2 \\&= \frac{x^2}{9 + 6 x} + 1 \quad\ge 1 \quad \text{since } x\ge -1. \end{align*}
The full Stirling approximation can similarly be derived from the series expansion: $$ z(1+1/x) = 1 + \frac{2z}{3} + \frac{z^2}{12} - \frac{2z^3}{135} + \frac{z^4}{864} + \dots $$ This form is nice as it allows easy integration with respect to the Gaussian pdf. E.g. $$ n\int_{-\infty}^{\infty}\frac{z^4}{864} e^{-n z^2/2} dz = \frac{\sqrt{2 \pi n}}{288n^{2}} $$ matching the expected $$ n! \sim \sqrt{2\pi n}\left(\frac{n}{e}\right)^n \left(1 +\frac{1}{12n}+\frac{1}{288n^2} - \frac{139}{51840n^3} -\frac{571}{2488320n^4}+ \cdots \right). $$
I think one can even give a uniform bound of $$ z(1+1/x) \le 1 + \frac{2z}{3} + \frac{z^2}{c}, $$ for some $c\approx 8$, strengthening the upper bound even more. However, the logarithmic inequalities become quite tedious.