Coefficients of the stirling's series expansion for the factorial. Knowing the Stirling's approximation for the Gamma function (factorial) for integers:
$$\Gamma(n+1)=n!\approx \sqrt{2\pi n}n^ne^{-n}\bigg(1+\frac{a_1}{n}+\frac{a_2}{n^2}+\cdots\bigg)$$
Using the above approximation one can write:
$$(n+1)!=\sqrt{2\pi(n+1)}(n+1)^{n+1}e^{-(n+1)}\bigg(1+\frac{a_1}{n+1}+\frac{a_2}{(n+1)^2}+\cdots\bigg)$$
We know that following recursion holds:
$$(n+1)!=(n+1)n!$$
One can rewrite this:
$$(n+1)!=(n+1)\sqrt{2\pi n}n^ne^{-n}\bigg(1+\frac{a_1}{n}+\frac{a_2}{n^2}+\cdots\bigg)$$
All this comes from: https://www.csie.ntu.edu.tw/~b89089/link/gammaFunction.pdf (Page 8-9)
Then the author gives this expansion to calculate the $a_k$ coefficients when $n$ becomes large.
Comparing these two expressions for $(n+1)!$ gives
$$1+\frac{a_1}{n}+\frac{a_2}{n^2}+\cdots=\bigg(1+\frac{1}{n}\bigg)^{n+1/2}e^{-1}\bigg(1+\frac{a_1}{n+1}+\frac{a_2}{(n+1)^2}+\cdots\bigg)$$
Then he says, that after "classical series expansion" this equals:
$$1+\frac{a_1}{n}+\frac{a_2-a_1+\frac{1}{12}}{n^2}+\frac{\frac{13}{12}a_1-2a_2+a_3+\frac{1}{12}}{n^3}+\cdots$$
I don't understand how he got there. Only thing that came to my mind was, that as $n\to\infty$ $\big(1+\frac{1}{n}\big)^n$ goes to $e$ but then we are left with $\big(1+\frac{1}{n}\big)^{1/2}$. When i expand this into binomial series, i get:
$$\bigg(1+\frac{1}{2n}-\frac{1}{8n^2}+\frac{1}{16n^3}\cdots\bigg)\cdot\bigg(1+\frac{a_1}{n+1}+\frac{a_2}{(n+1)^2}+\cdots\bigg)$$
And I'm stuck here. 
Or is there any other elementary way how to compute the coefficients $a_k$ of the stirling's series expansion for factorial/Gamma function?


We expand the series up to terms of $\frac{1}{n^3}$. We obtain for $\left|\frac{1}{n}\right|<1$
  \begin{align*}
\color{blue}{\left(1+\frac{1}{n}\right)^{n+\frac{1}{2}}e^{-1}}
&=e^{-1}e^{\left(n+\frac{1}{2}\right)\ln\left(1+\frac{1}{n}\right)}\\
&=\exp(-1)\exp\left[\left(n+\frac{1}{2}\right)\left(\frac{1}{n}-\frac{1}{2n^2}+\frac{1}{3n^3}+\cdots\right)\right]\\
&=\exp(-1)\exp\left(1+\frac{1}{12n^2}-\frac{1}{12n^3}+\cdots\right)\\
&=\exp\left(\frac{1}{12n^2}-\frac{1}{12n^3}+\cdots\right)\\
&\color{blue}{=1+\frac{1}{12n^2}-\frac{1}{12n^3}+\cdots}\tag{1}
\end{align*}

and applying the binomial series expansion

we obtain
  \begin{align*}
&\color{blue}{1+\frac{a_1}{n+1}+\frac{a_2}{(n+1)^2}+\frac{a_3}{(n+1)^3}+\cdots}\\
&\qquad =1+\frac{a_1}{n\left(1+\frac{1}{n}\right)}+\frac{a_2}{n^2\left(1+\frac{1}{n}\right)^2}
+\frac{a_3}{n^3\left(1+\frac{1}{n}\right)^3}+\cdots\\
&\qquad =1+\frac{a_1}{n}\left(1-\frac{1}{n}+\frac{1}{n^2}\cdots\right)
+\frac{a_2}{n^2}\left(1+\binom{-2}{1}\frac{1}{n}+\cdots\right)
+\frac{a_3}{n^3}\left(1+\cdots\right)\\
&\qquad\color{blue}{=1+a_1\cdot\frac{1}{n}+\left(-a_1+a_2\right)\frac{1}{n^2}
 +\left(a_1-2a_2+a_3\right)\frac{1}{n^3}+\cdots}\tag{2}\\ 
\end{align*}

Putting (1) and (2) together

we obtain
  \begin{align*}
&\left(1+\frac{1}{12n^2}-\frac{1}{12n^3}+\cdots\right)
\left(1+a_1\cdot\frac{1}{n}+\left(-a_1+a_2\right)\frac{1}{n^2}
 +\left(a_1-2a_2+a_3\right)\frac{1}{n^3}+\cdots\right)\\
&\quad\color{blue}{ =1+a_1\cdot\frac{1}{n}+\left(\frac{1}{12}-a_1+a_2\right)\frac{1}{n^2}
  +\left(-\frac{1}{12}+\frac{13}{12}a_1-2a_2+a_3\right)\frac{1}{n^3}+\cdots}
\end{align*}
and the claim follows.

