You can obtain the desired result using absolute convergence and some observations about the annulus of convergence of Laurent series.
If $\{x_n\}_{n \in \mathbb{Z}}$, then let $R_+(\{x_n\}) = ( \limsup_{n \to \infty} \sqrt[n]{|x_n|} )^{-1}$, and $R_-(\{x_n\}) = \limsup_{n \to \infty} \sqrt[n]{|x_{-n}|}$.
Two relevant results (all summations are over $\mathbb{Z}$):
(i) If $\sum_m \sum_n |a_{m,n}| < \infty$, then the summation can be rearranged. In particular, $\sum_m \sum_n a_{m,n} = \sum_k \sum_l a_{l,k-l}$. (Since $\phi(m,n) = (m,m+n)$ is a bijection of $\mathbb{Z}^2$ to $\mathbb{Z}^2$.)
(ii) If $\sum_{n} |x_n| < \infty$, then $R_+(\{x_n\}) \ge 1$. Similarly, $R_-(\{x_n\}) \le 1$.
Let $f(z) = \sum_n f_n (z-z_0)^n$, $g(z) = \sum_n g_n (z-z_0)^n$. Let $R_+ = \min(R_+(\{f_n\}),R_+(\{g_n\}))$, $R_- = \max(R_-(\{f_n\}),R_-(\{g_n\}))$.
Choose $R_-<r < R_+$. Then $\sum_m |f_m| r^m$ and
$\sum_n |g_n| r^n$ are absolutely convergent sequences, and so $\sum_m \sum_n |f_m||g_n| r^{m+n} < \infty$. From (i) we have $\sum_m \sum_n |f_m||g_n| r^{m+n} = \sum_k (\sum_l |f_l| |g_{k-l}|) r^k < \infty$.
If we let $c_k = \sum_l f_l g_{k-l}$, this gives $\sum_k |c_k| r^k < \infty$, and so
$R_+(\{c_kr^k\}) = \frac{1}{r} R_+(\{c_k \}) \ge 1$, or, in other words, $R_+(\{c_k \}) \ge r$. Since $r<R_+$ was arbitrary, it follows that $R_+(\{c_k \}) \ge R_+$. The same line of argument gives $R_-(\{c_k \}) \le R_-$.
Hence we have $c(z) = f(z)g(z) = \sum_n c_n (z-z_0)^k$ on $R_- < |z-z_0| < R_+$, where $c_k = \sum_l f_l g_{k-l}$.
For any given function, the Laurent series about a given point is unique. So you can pretty much use any method you want to find it. Using Cauchy integrals in the fashion you have described may work theoretically, but I'd imagine it would be really tedious to calculate.
The solution simply calculated the taylor series of $\sinh(z)$ and went from there.
Best Answer
Let's first consider $f(z)$ analytic inside and on a boundary $C$ of a simply connected region. Then Cauchy's integral formula is $f(z_0) = \frac{1}{2i\pi}\oint_C\frac{f(z)}{z-z_0}dz$.
Proof: Let $\Gamma$ be a circle of radius $\epsilon > 0$ inside $C$ with center $z_0$. Then $\oint_C\frac{f(z)}{z-z_0}dz = \oint_{\Gamma}\frac{f(z)}{z-z_0}dz$. The equation for a circle with center $z_0$ and radius $\epsilon$ is $\lvert z - z_0\rvert = \epsilon\iff z - z_0 = \epsilon e^{i\theta}$. Let $z = z_0 + \epsilon e^{i\theta}$ so $dz = i\epsilon e^{i\theta}d\theta$. \begin{align} \oint_{\Gamma}\frac{f(z)}{z-z_0}dz &= \int_0^{2\pi}\frac{f(z_0 + \epsilon e^{i\theta})}{\epsilon e^{i\theta}}i\epsilon e^{i\theta}d\theta\\ &= i\lim_{\epsilon\to 0}\int_0^{2\pi}f(z_0+\epsilon e^{i\theta})d\theta\\ &= 2i\pi f(z_0) \end{align} Thus, $f(z_0) = \frac{1}{2i\pi}\oint_C\frac{f(z)}{z-z_0}dz$.
Let $f(z)$ be analytic in a region bounded by two concentric circles $C_1$ and $C_2$ where the radius of $C_1$ is greater than the radius of $C_2$. By Cauchy's integral formula, we are able to show that $f(z_0) = \frac{1}{2i\pi}\oint_{C_1}\frac{f(z)}{z-z_0}dz-\frac{1}{2i\pi}\oint_{C_2}\frac{f(z)}{z-z_0}dz$.
Laurent's theorem: If $f(z)$ is analytic inside and on the boundary of an annular region bounded by two concentric circles centered at $z_0$ with radii $r_1$ and $r_2$, then for all $z$ in the annular region $$ f(z) = \sum_{n=0}^{\infty}a_n(z-z_0)^n+\sum_{n=1}^{\infty}\frac{a_{-n}}{(z-z_0)^n} $$ where the coefficients are defined as \begin{align} a_n &= \frac{1}{2i\pi}\oint_{C_1}\frac{f(w)}{(w-z)^{n+1}}dw\\ a_{-n} &= \frac{1}{2i\pi}\oint_{C_2}\frac{f(w)}{(w-z)^{-n+1}}dw \end{align}
Proof: By Cauchy's integral formula, we have $$ f(z) = \frac{1}{2i\pi}\oint_{C_1}\frac{f(w)}{w-z}dz-\frac{1}{2i\pi}\oint_{C_2}\frac{f(w)}{w-z}dz\tag{1} $$ for $z$ in the annular region. Since you are concerned with the $a_{-1}$, I am only going to show the work for the second integral in equation $(1)$. Consider $\frac{-1}{w-z}$. \begin{align} \frac{-1}{w-z} &= \frac{1}{z-z_0}\frac{1}{1 - \frac{w-z_0}{z-z_0}}\\ &= \frac{1}{z-z_0}\sum_{k=0}^{\infty}\Bigl(\frac{w-z_0}{z-z_0}\Bigr)^k\\ &= \frac{1}{z-z_0} + \frac{w-z_0}{(z-z_0)^2} + \frac{(w-z_0)^2}{(z-z_0)^3}+\cdots + \frac{(w-z_0)^{n-1}}{(z-z_0)^n}\\ &+ \text{higher order terms}\\ -\frac{1}{2i\pi}\oint_{C_2}\frac{f(w)}{w-z}dz &= \frac{1}{2i\pi}\oint_{C_2}\frac{f(w)}{z-z_0}dw + \frac{1}{2i\pi}\oint_{C_2}\frac{w-z_0}{(z-z_0)^2}f(w)dw + \frac{1}{2i\pi}\oint_{C_2}\frac{(w-z_0)^2}{(z-z_0)^3}f(w)dw\\ &+\cdots + \frac{1}{2i\pi}\oint_{C_2}\frac{(w-z_0)^{n-1}}{(z-z_0)^n}f(w)dw + \text{higher order terms}\\ &= \frac{a_{-1}}{z-z_0} + \frac{a_{-2}}{(z-z_0)^2} + \frac{a_{-3}}{(z-z_0)^3}+\cdots + \frac{a_{-n}}{(z-z_0)^n}\\ &+ \text{higher order terms} \end{align} Then $$ a_{-1}=\frac{1}{2i\pi}\oint_{C_2}f(w)dw. $$