I have just learned about the residue theorem so when I encountered the integral $$ I = \int_{0}^{+\infty} \frac{x^3}{e^x-1}dx $$ I tried to evaluate it with this method, which failed but produced something interesting.
Define the complex logarithm on $\mathbb C \setminus [0,+\infty)$ and consider the function $$ f(z) = \frac{z^3 \log z}{e^z-1} $$ with the path $\gamma$ given by this drawing
Provided that both $\oint_{C_\epsilon} fdz$ and $ \oint_{C_R}fdz$ goes to zero as $R\to\infty$ and $\epsilon \to 0$ the value of the original integral is given by the residue theorem
$$ I = -\frac{1}{2\pi i}\oint_\gamma fdz = -\sum_{k} \text{Res}(f, z_k) $$
$f$ has simple poles for $e^z=1$, $z=2k\pi i$ where $k \in \mathbb Z$ so the residues are
\begin{align*} \text{Res}(f, z_k) &= 0 && k=0 \\ \text{Res}(f, z_k) &= z^3 \log{z} \big|_{z=2k\pi i} = (2k\pi)^3 \log (2k\pi i) && k\neq 0 \end{align*}
The sum of the residue can be written as $$ \sum_{k=1}^{+\infty} (a_k + a_{(-k)}) $$
where \begin{align*} a_k &= i(2k\pi)^3\left(\log(2k\pi) + i\frac{\pi}{2}\right) \\ a_{(-k)} &= i(2k\pi)^3\left(-\log(2k\pi) - i\frac{3\pi}{2}\right) \end{align*}
Cancelling the terms out $$ I = \sum_{k=1}^{+\infty} i(2\pi k)^3(-i\pi) = 8\pi^4 \sum_{k=1}^{+\infty} k^3 $$
which diverges. Clearly something is wrong: maybe the integral on the circumference is not zero? I'm not sure. However if I replace the infinite sum with $\zeta(-3)=\frac{1}{120}$ I actually get the correct result $I = \frac{\pi^4}{15}$. Why does this work?
P.S. I managed to evaluate it in a more familiar way without involving complex analysis.
Let $f_s(z) = \frac {\exp(-s\log(z))\log (z)}{\exp(z)-1}$ and $\gamma$ be a keyhole curve going twice along the positive real axis and making a small loop around $0$, picking the branch of $\log$ where $\log(i) = i\pi/2$ (and so $\log(-i) = 3i\pi/2$ and not $-i\pi/2$ because the branch cut is on $\Bbb R^+$)
Since $f_s(z)$ gets smaller exponentially when $\Re(z)$ gets large, $J(s) = \int_\gamma f_s(z) dz$ is a well-defined function, and is in fact an entire function in $s$.
If $\Re(s) < 0$, then $f_s(z) \to 0$ when $z$ gets closer to $0$, so $J(s)$ can be evaluated by computing the difference between the two horizontal segments, and so $J(s) = \int_0^\infty \frac {\exp(-s\log(z))\log (z) - \exp(-s(\log(z)+2i\pi))(\log (z)+2i\pi)}{\exp(z)-1} dz$.
And so $J(-3) = \int_0^\infty \frac { -2i\pi \exp(3\log(z))}{\exp(z)-1} dz = -2i\pi \int_0^\infty \frac {z^3}{\exp(z)-1} dz$.
On the other hand, when $\Re(s) > 1$, we can approximate $J(s)$ by truncating the segment and closing with a large circle.
If $R_n = (2n+1)\pi$, then $\exp(R_n\exp(i\theta))-1$ stays away from $0$ independantly of $n$, and so $f_s(R_n\exp(i\theta)) = \exp(-s\log R_n-is\theta)(\log R_n +i \theta) / (\exp(z) - 1) = O(n^{- \Re(s)})$.
And so, the integral of $f_s dz$ on the large circle is a $O(n^{1-\Re(s)})$ and converges to $0$, and therefore
$J(s) = \lim _{n \to \infty} \int_{\gamma_n} f_s(z) dz$ where $\gamma_n$ is the closed curve you described with a radius $R_n$.
Now we can use the residue theorem to get $J(s) = 2i\pi \sum_1^\infty \exp(-s\log(2ik\pi))\log(2ik\pi) + \exp(-s(\log(2ik\pi)+i\pi))(\log(2ik\pi)+i\pi) \\ = 2i\pi \exp(-s\log(2i\pi))(\zeta'(s) + \zeta(s)\log(2i\pi)) + \\ 2i\pi \exp(-s(\log(2i\pi)+i\pi))(\zeta'(s) + \zeta(s)(\log(2i\pi)+i\pi))$
Now, we use uniqueness of analytic continuation to deduce that this equality is valid for all $s \neq 1$, and so
$J(-3) =\\ 2i\pi \exp(3\log(2i\pi))(\zeta'(-3) + \zeta(-3)\log(2i\pi)) + \\ 2i\pi \exp(3(\log(2i\pi)+i\pi))(\zeta'(-3) + \zeta(-3)(\log(2i\pi)+i\pi)) \\ = 2i\pi (\exp(3\log(2i\pi))(\zeta'(-3) + \zeta(-3)\log(2i\pi) - \zeta'(-3) - \zeta(-3)(\log(2i\pi)+i\pi))) \\ = - 2i\pi (\exp(3\log(2i\pi))( \zeta(-3)i\pi) \\ = -2i\pi (2i\pi)^3 i\pi \zeta(-3) = (-2i\pi)(8\pi^4 \zeta(-3)) $
Finally, putting both results together we get $\int_0^\infty \frac {z^3}{\exp(z)-1} dz = 8\pi^4 \zeta(-3)$