I understand how to find the residue for the simple pole at $z = 0$, but I don't understand how to find the residues at subsequent singularities located at $z = n\pi$. The solution tells me to set $w = z + n\pi$, which gives:
$$ z/\sin^2(z) = w/\sin^2(w) + n\pi/\sin^2(w) .$$
Which I understand how to get algebraically, but I don't understand where the conclusion that the residues of $z = n\pi$ are all equal to $1$ comes from. Is there also a way to do this by: $$ res(n\pi) = \lim_{z\to0}[(z - n\pi)(z/\sin^2(z))],$$ or is this not possible?
I study physics, so some simpler solutions would be helpful but any help is appreciated.
At $z=n\pi,n \neq 0$ you have a double pole, so the residue can be computed by the formula $\frac{1}{(n-1)!} \lim_{z \to n \pi} \frac{d^{n-1}}{dz^{n-1}} \left ( (z-n\pi)^n \frac{z}{\sin(z)^2} \right )$ with $n=2$, or in other words $\lim_{z \to n\pi} \frac{d}{dz} \left ( (z-n\pi)^2 \frac{z}{\sin(z)^2} \right )$.
This formula is rather standard, but if you don't see why it's different from the simple pole case, here's why it works. You start with $\sum_{k=-2}^\infty a_k (z-a)^k$. To get the residue at $a$, you want $a_{-1}$. Multiplying gets you to $\sum_{k=-2}^\infty a_k (z-a)^{k+2}$. Differentiating gets you to $\sum_{k=-2}^\infty (k+2) a_k (z-a)^{k+1}$. Now notice that the $k=-2$ term, which has a negative exponent, also has a coefficient of zero. As a result, the limit is just the coefficient of $(z-a)^0$, i.e. the coefficient in the $k=-1$ term, which is $a_{-1}$, exactly what you want.
Now why does the simple pole formula not work? This is because the $a_{-2}$ term would cause the limit that you would use in the simple pole case to not even exist.
With higher order poles, you have to shift $n$ times so that all the singular terms have nonnegative exponents, and then differentiate $n-1$ times so that the $n-1$ undesired singular terms all drop out, leaving just the one you want.