The Laurent series of the digamma function at the negative integers To find the Laurent series of  $\psi(z)$ at  $z= 0$, I would first find the Taylor series of $\psi(z+1)$ at $z=0$ and then use the functional equation of the digamma function.
Specifically,
$$\begin{align} \psi(z + 1) = \frac{1}{z} + \psi(z) &= \psi(1) + \psi'(1)z + \mathcal{O}(z^{2}) \\ &= -\gamma + \zeta(2) + \mathcal{O}(z^{2}) \\ &= -\gamma + \frac{\pi^{2}}{6}z + \mathcal{O}(z^{2}) \end{align}$$
which implies
$$ \psi(z) = -\frac{1}{z} - \gamma + \frac{\pi^{2}}{6} z + \mathcal{O}(z^{2}).$$
But I'm having trouble finding the Laurent series of $\psi(z)$ at the negative integers.
Since $\psi(z)$ has simple poles at the negative integers with residue $-1$, I know that the first term of the series must be $ \displaystyle\frac{-1}{z+n}$.
But I would like to determine more terms in the series.
EDIT: 
The series appears to be $$\begin{align} \psi(z) &= - \frac{1}{z+n} + (H_{n} - \gamma)+ \Big( H_{n}^{(2)} + \zeta(2) \Big) (z+n) + \Big( H_{n}^{(3)} - \zeta(3) \Big) (z+n)^{2} \\ &+  \Big( H_{n}^{(4)} + \zeta(4) \Big) (z+n)^{3}  + \Big( H_{n}^{(5)} - \zeta(5) \Big) (z+n)^{4} + \ldots \end{align}$$
SECOND EDIT: 
We can find the constant term by evaluating $\lim_{z \to -n} \Big( \psi(z) + \frac{1}{z+n} \Big)$.
Since $$ \begin{align} &\psi(z) + \frac{1}{z+n} \\ &= \psi(z+1) - \frac{1}{z} + \frac{1}{z+n} \\ &= \psi(z+2) - \frac{1}{z+1} - \frac{1}{z} + \frac{1}{z+n} \\ &= \ ... \ = \psi(z+n+1) - \frac{1}{z+n} - \frac{1}{z+n-1} - \ldots - \frac{1}{z+1} - \frac{1}{z} + \frac{1}{z+n}, \end{align}$$
we have $$ \lim_{z \to -n} \Big( \psi(z) + \frac{1}{z+n} \Big) = \psi(1) + 1 + \frac{1}{2} + \frac{1}{3} + \ldots + \frac{1}{n} = \psi(1) + H_{n} = H_{n} - \gamma.$$
Similarly, we can find the coefficient of the $(z+n)$ term by evaluating $ \lim_{z \to -n} \frac{\psi(z) + \frac{1}{z+n} - H_{n} + \gamma}{z+n}.$
Since $$ \psi_{1}(z) = \psi_{1}(z+n+1) + \frac{1}{(z+n)^{2}} + \frac{1}{(z+n-1)^{2}} +  \ldots + \frac{1}{z^2}, $$
we have $$ \lim_{z \to -n} \frac{\psi(z) + \frac{1}{z+n} - H_{n} + \gamma}{z+n} = \lim_{z \to -n} \Big(\psi_{1}(z) - \frac{1}{(z+n)^{2}} \Big) = \psi_{1}(1) + H_{n}^{(2)} =  H_{n}^{(2)} + \zeta(2) .$$
And we can find the coefficient of $(z+n)^{2}$ by evaluating $ \lim_{z \to -n} \frac{\psi(z) + \frac{1}{z+n} - H_{n} + \gamma -\big( H_{n}^{(2)} + \zeta(2) \big) (z+n)}{(z+n)^{2}} . $
Since $$\psi_{2}(z) = \psi_{2} (z+n+1) - \frac{2}{(z+n)^{3}} - \frac{2}{(z+n-1)^{3}} - \ldots - \frac{2}{z^{3}},$$
we have $\begin{align} \lim_{z \to -n} \frac{\psi(z) + \frac{1}{z+n} - H_{n} + \gamma -\big( H_{n}^{(2)} + \zeta(2) \big) (z+n)}{(z+n)^{2}} &= \lim_{z \to -n} \frac{\psi_{1}(x) - \frac{1}{(z+n)^{2}} -H_{n}^{(2)} - \zeta(2)}{2(z+n)} \\ &= \lim_{z \to -n} \frac{\psi_{2}(x) + \frac{2}{(z+n)^{3}}}{2} \\ &= \frac{1}{2} \Big( \psi_{2}(1) + 2 H_{n}^{(3)} \Big) \\ &= H_{n}^{(3)} - \zeta(3). \end{align}$
$ $
And so on.
 A: A few terms for starters....
$$
\Psi \left( x \right) = -{x}^{-1}-\gamma+1/6\,{\pi }^{2}x-\zeta 
 \left( 3 \right) {x}^{2}+{\frac {1}{90}}\,{\pi }^{4}{x}^{3}-\zeta 
 \left( 5 \right) {x}^{4}+{\frac {1}{945}}\,{\pi }^{6}{x}^{5}+O
 \left( {x}^{6} \right) 
$$
$$
\Psi \left( x \right) =- \left( x+1 \right) ^{-1}+1-\gamma+ \left( 1+
1/6\,{\pi }^{2} \right)  \left( x+1 \right) + \left( 1-\zeta  \left( 3
 \right)  \right)  \left( x+1 \right) ^{2}+ \left( {\frac {1}{90}}\,{
\pi }^{4}+1 \right)  \left( x+1 \right) ^{3}+ \left( 1-\zeta  \left( 5
 \right)  \right)  \left( x+1 \right) ^{4}+ \left( {\frac {1}{945}}\,{
\pi }^{6}+1 \right)  \left( x+1 \right) ^{5}+O \left(  \left( x+1
 \right) ^{6} \right) 
$$
$$
\Psi \left( x \right) =- \left( x+2 \right) ^{-1}+{\frac {3}{2}}-
\gamma+ \left( {\frac {5}{4}}+1/6\,{\pi }^{2} \right)  \left( x+2
 \right) + \left( {\frac {9}{8}}-\zeta  \left( 3 \right)  \right) 
 \left( x+2 \right) ^{2}+ \left( {\frac {1}{90}}\,{\pi }^{4}+{\frac {
17}{16}} \right)  \left( x+2 \right) ^{3}+ \left( {\frac {33}{32}}-
\zeta  \left( 5 \right)  \right)  \left( x+2 \right) ^{4}+ \left( {
\frac {1}{945}}\,{\pi }^{6}+{\frac {65}{64}} \right)  \left( x+2
 \right) ^{5}+O \left(  \left( x+2 \right) ^{6} \right) 
$$
$$
\Psi \left( x \right) =- \left( x+3 \right) ^{-1}+{\frac {11}{6}}-
\gamma+ \left( {\frac {49}{36}}+1/6\,{\pi }^{2} \right)  \left( x+3
 \right) + \left( {\frac {251}{216}}-\zeta  \left( 3 \right)  \right) 
 \left( x+3 \right) ^{2}+ \left( {\frac {1}{90}}\,{\pi }^{4}+{\frac {
1393}{1296}} \right)  \left( x+3 \right) ^{3}+ \left( {\frac {8051}{
7776}}-\zeta  \left( 5 \right)  \right)  \left( x+3 \right) ^{4}+
 \left( {\frac {1}{945}}\,{\pi }^{6}+{\frac {47449}{46656}} \right) 
 \left( x+3 \right) ^{5}+O \left(  \left( x+3 \right) ^{6} \right) 
$$
A: We will use the Laurent expansion around $x=0$
$$\psi_0(x)=-\frac{1}{x}-\gamma-\sum_{k\geq 1}\zeta(k+1)(-x)^k$$
Look the proof here 
Then
$$\tag{1}\psi_0(x+N)=-\frac{1}{x+N}-\gamma-\sum_{k\geq 1}(-1)^k\zeta(k+1)(x+N)^k$$
Now we use that 
$$\tag{2} \psi_0(x+N)=\psi_0(x)+\sum_{k=0}^{N-1}\frac{1}{x+k}$$
Now we look at the finite sum 
\begin{align}\sum_{k=0}^{N-1}\frac{1}{x+k}&=\sum_{k=0}^{N-1}\frac{1}{k-N+x+N}\\&=\sum_{k=0}^{N-1}\frac{1}{k-N}\frac{1}{1+\frac{x+N}{k-N}}\\ &= \sum_{k=0}^{N-1}\frac{1}{k-N}\sum_{m \geq 0}(-1)^m\left(\frac{x+N}{k-N} \right)^m \\&= \sum_{m \geq 0}(-1)^m(x+N)^m \sum_{k=0}^{N-1}\frac{1}{(k-N)^{m+1}}\\&=  -\sum_{m \geq 0}  H^{(m+1)}_{N} (x+N)^{m}\end{align}
Hence we have 
$$\psi_0(x+N)=\psi_0(x)-\sum_{m \geq 0} H^{(m+1)}_{N} (x+N)^{m}$$
Substitute in (1) to obtain 
$$\psi_0(x)-\sum_{m \geq 1}  H^{(m+1)}_{N} (x+N)^{m}=-\frac{1}{x+N}-\gamma-\sum_{k\geq 1}(-1)^k\zeta(k+1)(x+N)^k$$
$$\psi_0(x)=-\frac{1}{x+N}-\gamma-\sum_{k\geq 1}(-1)^k\zeta(k+1)(x+N)^k+\sum_{m \geq 0} H^{(m+1)}_{N} (x+N)^{m}$$
$$\psi_0(-x)=\frac{1}{x-N}-\gamma-\sum_{k\geq 1}\zeta(k+1)(x-N)^k+\sum_{m \geq 0}(-1)^m H^{(m+1)}_{N} (x-N)^{m}$$
$$\psi_0(-x)+\gamma=\frac{1}{x-N}+H_N+\sum_{k\geq 1}((-1)^kH^{(k+1)}_{N}-\zeta(k+1))(x-N)^k$$

Proof of (2)
$$\psi_0(x+N)=\psi_0(x)+\sum_{k=0}^{N-1}\frac{1}{x+k}$$
By induction for $N=1$ we get 
$$\psi_0(x+1)=\psi_0(x)+\frac{1}{x}$$
Which is true and can be proved using 
$$\psi_0(x+1)=-\gamma+\sum_{n\geq 1}\frac{x}{n(n+x)}$$
Now for the inductive step , assume 
$$\psi_0(x+N)=\psi_0(x)+\sum_{k=0}^{N-1}\frac{1}{x+k}$$
Then 
$$\psi_0(x+N+1)=\psi_0(x+N)+\frac{1}{x+N+1}=\psi_0(x)+\sum_{k=0}^{N}\frac{1}{x+k}$$
