Question: I want to prove, the Laurent expansion of gamma function.
\begin{align}
\Gamma(z) = \frac1z-\gamma+\frac12\left(\gamma^2+\frac {\pi^2}6\right)z-\frac16\left(\gamma^3+\frac {\gamma\pi^2}2+2 \zeta(3)\right)z^2+O(z^3).
\end{align}

First, my guess of obtaing above expansion, 
is starting from the definitions of gamma function
\begin{align}
\Gamma(z) &= \int\_0^{\infty} dt e^{-t} t^{z-1} \\
& = \int\_1^\infty dt e^{-t}t^{z-1} + \int\_0^1dt e^{-t} t^{z-1} \\
& = \int\_1^\infty dte^{-t}t^{z-1} + \int\_0^1 dt t^{z-1} \sum\_{n=0}^{\infty}\frac{(-1)^n}{n!}t^n \\
& = \int\_1^\infty dt e^{-t}t^{z-1} + \sum\_{n=0}^{\infty} \frac{(-1)^n}{n!}\frac{1}{z+n}
\end{align}
This only gives the gamma function as a function of $\frac{1}{z}$...

Or should I start with 
\begin{align}
\Gamma(z) = \lim\_{n \rightarrow \infty} \frac{n! n^z}{z(z+1) \cdots(z+n)}
\end{align}

Answer: 
Let $\Gamma(z)$ be represented by the integral

$$\Gamma(z)=\int\_0^\infty x^{z-1}e^{-x}\,dx\tag1$$

for $\text{Re}(z)>0$. Integrate by parts the integral in $(1)$ with $u=e^{-x}$ and $v=\frac{x^z}{z}$ reveals

$$\Gamma(z)=\frac1z\int\_0^\infty x^ze^{-x}\,dx\tag2$$

Next, we expand $x^z$ in a power series of $z$ to obtain

$$\begin{align}
\Gamma(z)&=\frac1z\sum\_{n=0}^\infty \frac{z^n}{n!}\int\_0^\infty e^{-x}\log^n(x)\,dx\\\\
&=\frac1z+\underbrace{\int\_0^\infty \log(x)e^{-x}\,dx}\_{=-\gamma}+\frac12 z\underbrace{\int\_0^\infty \log^2(x)e^{-x}\,dx}\_{\gamma^2+\frac{\pi^2}6}+\frac16z^2 \underbrace{\int\_0^\infty \log^3(x)e^{-x}\,dx}\_{-\gamma^3-\gamma\pi^2/2-2\zeta(3)}+O(z^3)
\end{align}$$

as was to be shown.

---

**NOTE:**

The coefficients $\int\_0^\infty \log^n(x)e^{-x}\,dx$ for $n=2,3$ can be found by using the relationship, $\Gamma'(x)=\Gamma(x)\psi(x)$, between the logarithmic derivative of the Gamma function and the Digamma function, along with values of $\psi'(1)=\zeta(2)$ and $\psi''(1)=-2\zeta(3)$.
An idea for you to develop:

The Weierstrass Formula tells us that

$$\frac1{\Gamma(z)}=ze^{\gamma z}\prod\_{n=1}^\infty\left(1+\frac zn\right)e^{-z/n}$$

Now take logarithms on both sides to get a more or less well known relation:

$$-\log\Gamma(z)=\log z+\gamma z+\sum\_{n=1}^\infty\left[\log\left(1+\frac zn\right)-\frac zn\right]$$

Now differentiate the above to get the logarithmic derivative of the Gamma function:

$$\frac{\Gamma'(z)}{\Gamma(z)}=-\frac1z-\gamma-\sum\_{n=1}^\infty\frac1n\left[\frac n{z+n}-1\right]=-\frac1z-\gamma+\sum\_{n=1}^\infty\frac z{n(z+n)}$$

and etc. You can try to integrate the $\;-\dfrac1z\;$ term into the series, too.