"In general, when we have more than a single pole, i.e a double pole, triple pole etc. is it best to just try and find the Laurent expansion about the pole?"
In my experience, for small orders, usually it's easier to use a well known formula.
Given an holomorphic function $\varphi$ with a pole order $k$ at $z_0$, one has, around $z_0$,
$$\varphi (z)=\sum \limits_{n=-k}^\infty \left(c_n(z-z_0)^n\right),$$ for a certain complex sequence $\left(c_n\right)_{n\in \mathbb Z}$, where $\forall n\in \mathbb Z(n<-k\implies c_n=0)$.
Multiply by $(z-z_0)^k$ to get
$$\begin{align}(z-z_0)^k\varphi(z)&=\sum \limits_{n=-k}^\infty \left(c_n(z-z_0)^{n+k}\right)\\ &=c_{-k}+c_{-k+1}(z-z_0)+\ldots +\underbrace{c_{-1}}_{\text{Res}(\varphi , z_0)}(z-z_0)^{-1+k}+\ldots .\end{align}$$
Differentiating $k-1$ times yields $$\dfrac{\mathrm d^{k-1}}{\mathrm dz^{k-1}}\left(z\mapsto (z-z_0)^k\varphi(z)\right)(z)=(k-1)!c_{-1}+k!c_0(z-z_0)+\ldots.$$
Thus, evaluating the above at $z_0$ gives you the residue since everything on the RHS except for $(k-1)!c_{-1}$ goes astray. However $\varphi$ is not defined at $z_0$, so you can't really evaluate at $z_0$, that's OK, take the limit at $z_0$ which exists by definition of pole of order $k$.
This justifies the formula $$\text{Res}(\varphi, z_0)=\dfrac 1{(k-1)!}\lim \limits_{z\to z_0}\left[\dfrac{\mathrm d^{k-1}}{\mathrm dz^{k-1}}\left(z\mapsto (z-z_0)^k\varphi(z)\right)(z)\right].$$
For this particular example, I think it's better to compute $c_{-1}$ by using the series expansion.
Set $g\colon \mathbb C\to \mathbb C, z\mapsto \cos(z)-1, S=\{2k\pi i\colon k\in \mathbb Z\}, h\colon \mathbb C\setminus S\to \mathbb C, z\mapsto \dfrac1{\left(e^z-1\right)^2}$, and $f\colon \mathbb C\setminus S\to \mathbb C, z\mapsto g(z)h(z)$.
Let $k\in \mathbb Z$ be such that $k\neq0$ and write $s=2k\pi i$.
'Around' $s$ each of the functions above is (holomorphic, hence) analytic.
Therefore, for some sequences $(a_n)_{n\in\mathbb Z}, (b_n)_{n\in \mathbb Z}, (c_n)_{n\in \mathbb Z}$ and for all $z$ in a sufficiently small punctured disc centered at $s$ one has
$$g(z)=\sum \limits_{n\in \mathbb Z}\left(a_n(z-s)^n\right),$$
$$h(z)=\sum \limits_{n\in \mathbb Z}\left(b_n(z-s)^n\right),$$
$$f(z)=\sum \limits_{n\in \mathbb Z}\left(c_n(z-s)^n\right).$$
You've determined that $s$ is a pole of order $2$ of $f$. The same happens to $h$ since $\lim \limits_{z\to s}\left(\left|(z-s)h(z)\right|\right)=\infty$ and $$\lim \limits_{z\to s}\left((z-s)^2h(z)\right)=\lim \limits_{z\to s}\left(\dfrac{2(z-s)}{2(e^z-1)e^z}\right)=\lim \limits_{z\to s}\left(\dfrac{1}{2e^{2z}-e^z}\right)=1\neq 0 \tag{$\huge ☺$}$$
Thus $$\forall n\in \mathbb Z(n<0\implies a_n=0),$$
$$\forall n\in \mathbb Z(n<2\implies b_n=0),$$
$$\forall n\in \mathbb Z(n<2\implies c_n=0).$$
Therefore $$\sum \limits_{n=-2}^\infty\left(c_n(z-s)^n\right)=f(z)=g(z)h(z)=\sum \limits_{n=0}^\infty\left(a_n(z-s)^n\right)\sum \limits_{n=-2}^\infty\left(b_n(z-s)^n\right).$$
Recall that the goal is to find $c_{-1}$.
In particular and more pictorially one has
$$\begin{align}
g(z)h(z)&=\left(a_0+a_1(z-s)+\ldots\right)\left(b_{-2}(z-s)^{-2}+b_{-1}(z-s)^{-1}+\ldots\right)\\
&=\left(a_0b_{-2}(z-s)^{-2}+(a_0b_{-1}+a_1b_{-2})(z-s)^{-1}+\ldots\right).
\end{align}$$
Consequently $c_{-1}=a_0b_{-1}+a_1b_{-2}$.
Because $g$ is analytic, $a_0$ and $a_1$ are just the coefficients of the taylor series of $g$ around $s$, that is, $$a_0=f(s)=\dfrac{e^{i\cdot 2k\pi i}+e^{-2k\pi i\cdot i}}{2}-1=\cosh(2k\pi)-1$$ and $$a_1=f'(s)=-\sin(s)=-\dfrac{e^{2k\pi i\cdot i}-e^{-2k\pi i\cdot i}}{2i}=\dfrac{e^{2k\pi}-e^{-2k\pi}}{2i}=-i\sinh(2k\pi).$$
Where as $\displaystyle b_{-2}=\dfrac{1}{2\pi i}\int _\gamma \dfrac{h(z)}{(z-s)^{-1}}\mathrm dz$ and $\displaystyle b_{-1}=\dfrac{1}{2\pi i}\int _\gamma h(z)\,\mathrm dz$, where $\gamma$ is an appropriate positively parametrized circle, centered at $s$, going around exactly once.
The integrals can be computed using the residue theorem.
By choosing $\gamma$ with a small enough radius, the only singularity inside the $\text{im}(\gamma)$ is $s$.
One gets $$\int _\gamma \dfrac{z-s}{(e^z-1)^2}\mathrm dz=2\pi i\text{Res}\left(z\mapsto \dfrac{z-s}{(e^z-1)^2}, s\right).$$
The formula given in the first part of this answer, together with $\huge ☺$, yields $\text{Res}\left(z\mapsto \dfrac{z-s}{(e^z-1)^2}, s\right)=1$, thus $b_{-2}=1$.
For the other, since $s$ is a pole of order $2$, one gets
$$\begin{align}
\int _\gamma h(z)\,\mathrm dz=\int _\gamma \dfrac{1}{(e^z-1)^2}\mathrm dz&=2\pi i\text{Res}(h,s)\\
&=2\pi i\lim \limits_{z\to s}\left[\dfrac{\mathrm d}{\mathrm dz}\left(z\mapsto (z-s)^2h(z)\right)(z)\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{\mathrm d}{\mathrm dz}\left(w\mapsto \dfrac{w^2}{(e^{w+2k\pi i}-1)^2}\right)(z)\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{\mathrm d}{\mathrm dz}\left(w\mapsto \dfrac{w^2}{(e^{w}-1)^2}\right)(w)\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{2w(e^w-1)^2-2w^2(e^w-1)e^w}{(e^w-1)^4}\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{2w(e^{w}-1)-2w^2e^w}{(e^{w}-1)^3}\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{2w\left(w+\frac{w^2}2+\ldots\right)-2w^2(1+w+\ldots)}{(w+\ldots)^3}\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{2w^2+w^3-2w^2-2w^3+\ldots}{w^3+\ldots}\right]\\
&=2\pi i\lim \limits_{w\to 0}\left[\dfrac{-w^3+\ldots}{w^3+\ldots}\right]\\
&=-2\pi i,
\end{align}$$
therefore $b_{-1}=-1$.
Thus $\text{Res}(f,s)=c_{-1}=1-\cosh(2k\pi)-i\sinh(2k\pi)$. This agrees with wolfram alpha.
Best Answer
Alright let's see if I can answer this. We have the function $$f(z)=\frac{z}{1-\cos(z)}$$ We want to find all poles and residues at said poles.
First look at the numerator. This is just the identity function, and will not have a pole as it is analytical throughout the complex plane. However, we can note that at the origin, it is $0$, which may cause problems. Let's set this aside for now.
For the denominator, the poles exist at the solutions of $$1-\cos(z)=0$$ $$\implies z=2\pi n,\quad n\in\mathbb{Z}$$
Now what about the degree of the poles?
We could observe that by reducing the entire denominator into one single term, to wit, $1-\cos(z)=2\sin^2\left(\frac z2\right)$, this would mean that the multiplicity of all poles should be two. Now at $0$, we can see that this would become an indeterminate form. You can intuitively think (similar to the $\frac{z}{\sin(z)}$ case) that the zero from the top cancels out one degree of the pole and makes the pole at the origin degree $1$ instead. I don't recommend doing this because for trigonometric functions the multiplicity of the pole based purely on power is a tentative thing.
What you did is also plausible. Recall that if $z=c$ is a pole of order $n$, then the residue of $f$ around $c$ can be found by $$\mathop{\mathrm{Res}}_{z = c}f(z) = \frac{1}{(n-1)!} \lim_{z \to c} \frac{\text{d}^{n-1}}{\text{d}z^{n-1}} \left[ (z-c)^n f(z) \right]$$ We can test degree by degree upwards, letting $n=1$ first, and increasing $n$ based on the following fact: if the limit is infinite, increase $n$ by $1$ and try again. By this, we test that all poles have a multiplicity of $2$ except the one at the origin which has a multiplicity of $1$. You can also find all residues this way(though I do indeed get 2, contrary to what you get but w/e)
This is of course tedious, so we can instead use a third method. Invert $f$ by raising it to the $-1$st power. Call this $g(z)$.
If $g(z) = g '(z) = 0$, $g ''(z)\ne0$ , then $z$ is a double zero.
If $g(z) = g '(z) = g ''(z) = 0$, $g '''(z)\ne0$ , then $z$ is a triple zero.
And so on.
By this, we can see that all zeroes have a multiplicity of $2$, except for the one at $0$ which has a multiplicity of $1$.
so the answer to this question (as already established in the comments is W|A is wrong and you're correct)
Alternatively, you could also use a laurent series expansion. We have $$1-\cos(z)=1-\left[\cos(2 n \pi) - \sin(2 n \pi) (z - 2n\pi) - \frac 12 \cos(2 n\pi) (z - 2 n\pi)^2 + \frac 16 \sin(2 n\pi)(z - 2 n\pi)^3+\frac 1{24}\cos(2n\pi)(z-2 n\pi)^4-\frac1{120}\sin(2n\pi)(z-2n\pi)^5-\frac1{720}\cos(2n\pi)(z-2n\pi)^6+\mathcal{O}(z-2n\pi)^7 \right]$$
so $$\frac{z}{1-\cos(z)}=\frac{z}{1-\cos(2 n \pi) + \sin(2 n \pi) (z - 2n\pi) + \frac 12 \cos(2 n\pi) (z - 2 n\pi)^2 - \frac 16 \sin(2 n\pi)(z - 2 n\pi)^3-\frac 1{24}\cos(2n\pi)(z-2 n\pi)^4+\frac1{120}\sin(2n\pi)(z-2n\pi)^5+\frac1{720}\cos(2n\pi)(z-2n\pi)^6+\mathcal{O}(z-2n\pi)^7}$$ There are two cases. First, $n\neq 0$. WLOG $n=1$. Plug this in and after some manipulation we see that the lowest two terms are $$\frac{4\pi}{(z-2\pi)^2} +\frac{2}{z-2\pi}+...$$
We can tell two things from this. First, the lowest degree is the multiplicity of the pole (which is multiplicity $2$). Next, the residue is the coefficient of the term with degree $-1$, and it is $2$.
The next case is $n=0$, to which we can plug this in, and a nice symmetry of the denominator and numerator will cause the series to become $$\frac2z+\frac z6+...$$ Hence, the degree of the pole is $1$, and the residue at this pole is still $2$.
As a matter of fact, if you do some tedious algebra bashing of the more general series, you will see that the coefficient of the term with the degree $-1$ will always be $2$(and hence, residue will always be $2$).
Usually when using laurent series however, we just focus on finding the residue, not the degree of the pole unless explicitly asked. Furthermore, as you can see in this case, it may be quite inconvenient, (and if you work out the algebra i've ommitted and shortcut mathematica used) especially if it is expanded in the denominator.
Typically building a Laurent series, we would use a geometric series, binomial series, taylor series(like in this case), or any easily availible series to expand out a known term(s) and simplify enough such that you can find the coefficient of the $z^{-1}$ term(ie matching indicies, shortcut PFD, etc).
the answer to this question would then be that if you don't want to use the higher order pole residue formula, find the function's laurent series and then look for the coefficient of the term $z^{-1}$