The way I've usually shown Stirling's Asymptotic Approximation starts with the formula \begin{align} n! &=\int_0^\infty x^ne^{-x}\;\mathrm{d}x\\ &=\int_0^\infty e^{-x+n\log(x)}\;\mathrm{d}x\\ &=\int_0^\infty e^{-nx+n\log(nx)}n\;\mathrm{d}x\\ &=n^{n+1}\int_0^\infty e^{-nx+n\log(x)}\;\mathrm{d}x\\ &=n^{n+1}e^{-n}\int_{-1}^\infty e^{-nx+n\log(1+x)}\;\mathrm{d}x\\ &=n^{n+1}e^{-n}\int_{-\infty}^\infty e^{-nu^2/2}x'\;\mathrm{d}u\tag{1} \end{align} Where $$u^2/2=x-\log(1+x)=x^2/2-x^3/3+\dots\quad$$ Thus, $$x$$ has a power series in $$u$$ like $$x=0+u+\dots\quad$$ Let's compute some of that series.
Note that $$u(1+x)=xx'\tag{2}$$ and suppose that $$x=\sum_{k=1}^\infty a_ku^k\tag{3}$$ where $$a_1=1$$. Then, looking at the coefficient of $$u^n$$ in $$(2)$$, we get \begin{align} a_{n-1} &=\sum_{k=1}^na_{n-k+1}ka_k\\ &=\sum_{k=1}^na_k(n-k+1)a_{n-k+1}\\ &=\frac{n+1}{2}\sum_{k=1}^na_ka_{n-k+1}\tag{4} \end{align} Therefore, we get the recursion $$a_n=\frac{a_{n-1}}{n+1}-\frac{1}{2}\sum_{k=2}^{n-1}a_ka_{n-k+1}\tag{5}$$ Thus, we get $$x=u+\tfrac{1}{3}u^2+\tfrac{1}{36}u^3-\tfrac{1}{270}u^4+\tfrac{1}{4320}u^5+\tfrac{1}{17010}u^6-\tfrac{139}{5443200}u^7+\dots\tag{6}$$ In $$(1)$$, we integrate $$x'$$ against an even function of $$u$$, so we only need the even order terms of $$x'$$. Thus, using the identity $$\int_{-\infty}^\infty e^{-u^2}u^{2n}\;\mathrm{d}u=\tfrac{(2n-1)!!}{2^n}\sqrt{\pi}\tag{7}$$ we get \begin{align} n! &=n^{n+1}e^{-n}\int_{-\infty}^\infty e^{-nu^2/2}\left(1+\tfrac{1}{12}u^2+\tfrac{1}{864}u^4-\tfrac{139}{777600}u^6+\dots\right)\;\mathrm{d}u\\ &=\sqrt{2n}\;n^ne^{-n}\int_{-\infty}^\infty e^{-u^2}\left(1+\tfrac{1}{6n}u^2+\tfrac{1}{216n^2}u^4-\tfrac{139}{97200n^3}u^6+\dots\right)\;\mathrm{d}u\\ &=\sqrt{2\pi n}\;n^ne^{-n}\left(1+\frac{1}{12n}+\frac{1}{288n^2}-\frac{139}{51840n^3}+\dots\right)\tag{8} \end{align}