Laurent series of $f(z) = \frac{z^2+1}{\sin(z)}$ How does one find the first four terms of the Laurent series of $f(z) = \frac{z^2+1}{\sin(z)}$? My approach was this:
$(z^2+1) = f(z)\sin(z) = \left(\sum\limits_{n=-\infty}^\infty c_n z^n \right)\sin(z)=\left(\sum\limits_{n=-\infty}^\infty c_n z^n \right)\left(z-\frac{z^3}{6}+\frac{z^5}{120}-\frac{z^7}{5040} \right)$.
How do we proceed further?
 A: For example, very close to zero we have
$$\frac{z^2+1}{\sin z}=\frac{z^2+1}{z-\frac{z^3}{3!}+\frac{z^5}{5!}-\ldots}=\frac{z^2+1}{z\left(1-\frac{z^2}6+\frac{z^4}{120}-\ldots\right)}\stackrel{\text{Devel. of geom. series}}=$$
$$=\left(z+\frac1z\right)\left(1+\frac{z^2}6+\frac{z^4}{36}+\ldots\right)=\frac1z+\frac76z+\frac7{36}z^3+\ldots$$
With the above you already have the pole at $\;z=0\;$ is a simple one and its residue is $\;1\;$ .
A: $$\frac{z^2+1}{\sin z}=\biggl( z+\frac 1z\biggr)\frac z{\sin z},$$
so all you have to find is the expansion of
$$\frac z{\sin z}=\frac 1{1-\cfrac{z^2}6+\cfrac{z^4}{120}+\cfrac{z^6}{5040}+o(z^7)}$$
It can be obtained with a division by increasing powers of $x$:
$$\begin{array}{rrrrr}
&&1&{}+\dfrac{z^2}6&{}+\dfrac{7z^4}{360}&{}+\dfrac{31z^6}{15120}\\
1-\dfrac{z^2}6+\dfrac{z^4}{120}-\dfrac{z^6}{5040}&\Big(&1\\%
&&-1&{}+\dfrac{z^2}6&{}-\dfrac{z^4}{120}&{}+\dfrac{z^6}{5040}\\
\hline 
&&&\dfrac{z^2}6&{}-\dfrac{z^4}{120}&{}+\dfrac{z^6}{5040}\\
&&&-\dfrac{z^2}6&{}+\dfrac{z^4}{36}&{}-\dfrac{z^6}{360}\\
\hline 
%&&&&{}+\dfrac{7z^4}{360}&{}-\dfrac{z^6}{840}\\
%&&&&{}-\dfrac{7z^4}{360}&{}\dfrac{7z^6}{2160}\\
%\hline
&&&&&\dfrac{31z^6}{15120}\end{array}$$ 
whence the Laurent series after multiplication by $z+\dfrac1z$:
\begin{alignat*}{6}
\frac{z^2+1}{\sin z}=\biggl(z+\frac 1z\biggr)\frac z{\sin z}&=&&z+\dfrac{z^3}6&&+\dfrac{7z^5}{360}+\dfrac{31z^7}{15120}+o(z^8)
\\
&=\quad \frac1z+&&\dfrac{z}6+\dfrac{7z^3}{360}&&+\dfrac{31z^5}{15120}+o(z^6)\\
& =
\quad \color{red}{\frac1z+}&&\color{red}{\dfrac{7z}6+\dfrac{67z^3}{360}}&&\color{red}{{}+\dfrac{65z^5}{3024}+o(z^6)}.\end{alignat*}
Added:  Division by increasing power order
It's like Euclidean division of polynomials $A(x)$ (dividend) and $B(x)$ (divisor), except it is defined when $B(0)\neq 0$, i.e. the divisor has a non-zero constant term, and at each step, one divides the lowest term of the dividend by the constant term of the divisor in order to eliminate the lowest term of the dividend.
The  result we use is the following:

Let $A(x)$, $B(x)$ be polynomials. For each natural number $n$, there exists a unique pair of polynomials $Q_n(x)$ and $R_n(x)$ such that
  $$A(x)=Q_n(x)B(x)+x^{n+1}R_n(x),\qquad\deg Q_n(x)\le n.$$

From unicity, we see that, for $m\le n$, $Q_m(x)$ is just the polynomial $Q_n(x)$, truncated at degree $m$.
This division can be used for the expansion in power series of the quotient $f(x)/g(x)$ of two analytic functions, and to find quickly the decomposition into partial fractions in the case of multiple poles.
A: Your approach is correct, but without determining the order of the pole $z=0$ you are stuck because each coefficient of $f \cdot \sin$ will, in principle, be a sum with infinitely many terms. Finding the order of the pole of $f$ will allow us to kill infinitely many of these, leaving only finitely many of them, as you are going to see.
Remember that the Taylor series of $\sin$ around $z=0$ begins with $z^1$, so $0$ is a zero of order $1$ for $\sin$, therefore it is a pole of order $1$ for $\frac 1 \sin$.
For even greater clarity, remember that $z_0$ is a pole of order $n>0$ for $f$ if and only if $n$ is smallest with the property $\lim \limits _{z \to z_0} (z-z_0)^n f(z) \ne 0$. Let's check this for your $f$ and $n=1$:
$$\lim _{z \to 0} z \ \frac {z^2 + 1} {\sin z} = \lim _{z \to 0} \frac z {\sin z} (z^2 + 1) = \lim _{z \to 0} \frac z {\sin z} \lim _{z \to 0} (z^2 + 1) = 1 \cdot 1 = 1 .$$
Is $n=1$ the smallest such strictly positive number? Yes, because below it there is only $0$, and $0$ is not allowed to be the order of a pole, by definition.
Having found the order of the pole of $f$ is essential, since this allows us to restrict the lower bound for the summation index:
$$f(z) = \sum _{n = -1} ^\infty c_n z^n$$
(notice that $n$ starts now at $-1$, not at $-\infty$).
Plugging this into your own approach (and keeping only terms up to and including $z^3$, because you want only the first 4 terms) gives
$$z^2 + 1 = \left( \sum _{n =-1} ^\infty c_n z^n \right) \left( z - \frac {z^3} {3!} + \frac {z^5} {5!} - \frac {z^7} {7!} + \dots \right) = \\
\left( c_{-1} - \frac {c_{-1}} {3!} z^2 + \dots \right) + \left( c_0 z - \frac {c_0} {3!} z^3 + \dots \right) + \left( c_1 z^2 + \dots \right) + \left( c_2 z^3 + \dots \right) + \dots ,$$
so equating the coefficients of the corresponding powers of $z$ in both hands of the equation gives
$$c_{-1} = 1, \quad c_0 = 0, \quad -\frac {c_{-1}} 6 + c_1 = 1, \quad -\frac {c_0} 6 + c_2 = 0 ,$$
whence it follows that
$$c_{-1} = 1, \quad c_0 = 0, \quad c_1 = \frac 7 6, \quad c_2 = 0 ,$$
so
$$f(z) = \frac 1 z + \frac 7 6 z + \dots .$$
Notice that there is no surprise in the fact that $c_0 = c_2 = 0$, because
$$f(-z) = \frac {(-z)^2 + 1} {\sin (-z)} = \frac {z^2 + 1} {-\sin z} = -f(z) ,$$
so $f$ is odd, and the coefficients of even order of an odd function are always $0$.
