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?

$$\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.
