Contour integral for $x^3/(e^x-1)$?

I assume by "or is this going to need some something else" you meant that you are open to non-contour integration techinques, I will show you one of those and let the contour integration be left to someone else (it's not too bad, I just wish I had graphing software).

The classic way to evaluate this integral is as follows

$$\begin{aligned}\int_0^\infty\frac{x^m}{e^x-1}\;dx &= \int_0^{\infty}\frac{e^{-x}x^m}{1-e^{-x}}\;dx\\ &=\int_0^{\infty}\sum_{n=0}^{\infty}x^me^{-x(n+1)}\;\\ &= \sum_{n=0}^{\infty}\int_0^{\infty}x^me^{-x(n+1)}\;dx\\ &=\Gamma(m+1)\sum_{n=0}^{\infty}\frac{1}{(n+1)^{m+1}}\\ &=\Gamma(m+1)\zeta(m+1)\end{aligned}$$

This is, in fact, the way one defines $\zeta(x)$ for $x>1$.

To solve the integral I thought of a rectangle of vertices $0$, $R$, $R+i2\pi $ and $i2\pi$. Since the function we are going to use is singular in $i2k\pi $ we have to indent the first and the last vertex by using a quarter of circle for each one: let's make them small, of radius $\epsilon$.

Calling this contour $\Gamma$ we have:

$$ \oint_{\Gamma} \frac{z^4}{e^z -1}\mathrm{d}x = 0 $$ by Cauchy's theorem.

As you will notice we have by now 4 segments and 2 quarters of circle to integrate on, namely: $$ \int_\epsilon ^R \frac{x^4}{e^x - 1}\mathrm{d}x + \int_0 ^{2\pi} \frac{(R+iy)^4}{e^{R+iy}- 1}i\mathrm{d}y + \int_{R}^\epsilon \frac{(x+i2\pi)^4}{e^{x+i2\pi}-1}\mathrm{d}x + \int_0 ^{-\frac{\pi}{2}}\frac{(2 \pi i + \epsilon e^{i\theta})^4}{e^{2\pi i + \epsilon e^{i\theta}}-1} i\epsilon e^{i\theta}\mathrm{d}\theta + \int_{2\pi- \epsilon}^{\epsilon}\frac{(iy)^4}{e^{iy}-1}i\mathrm{d}y+ \int_{\frac{\pi}{2}}^{0}\frac{(\epsilon e^{i\theta})^4}{e^{\epsilon e^{i\theta}}-1}i\epsilon e^{i\theta} \mathrm{d}\theta = 0. $$

We can first expand the power in the third integral and note that its first term cancels out with the first integral.

Before writing anything else we take into account the limits $R \to +\infty$, which cancels out the second integral, and $\epsilon \to 0^+$, which annuls the last one and yields a finite value $-8i\pi^5$ in the 4th one.

Therefore we get $$ -i8\pi \int_0 ^\infty \frac{x^3}{e^x - 1} \mathrm{d}x + 24\pi^2\int_0 ^\infty \frac{x^2}{e^x -1}\mathrm{d}x + i 32 \pi^3 \int_0 ^\infty \frac{x}{e^x - 1}\mathrm{d}x- 16\pi^4\int_0 ^\infty \frac{1}{e^x - 1}\mathrm{d}x -i8\pi^5+\frac{i}{2} \int_0 ^{2\pi} y^4 \mathrm{d}y - \frac{1}{2} \int_0 ^{2\pi} \frac{y^4 \sin y}{1-\cos y}\mathrm{d}y=0. $$ Where we have split the 5th integral into its real and imaginary parts.

By considering the expression of the imaginary parts only: $$ -8\pi \int_0 ^\infty \frac{x^3}{e^x - 1} \mathrm{d}x + 32\pi^3 \int_0 ^\infty \frac{x}{e^x - 1}\mathrm{d}x - 8\pi^5 + \frac{16}{5}\pi^5 = 0. $$

The second integral gives $\frac{\pi^2}{6}$ (for a solution of this one a contour similar to the one we've used here is needed; I'm pretty sure it has already been solved here on Math.SE).

Finally: $$ \int_0 ^\infty \frac{x^3}{e^x - 1} \mathrm{d}x = \frac{\pi^4}{8} \left(\frac{16}{3}-8+\frac{16}{5}\right) = \frac{\pi^4}{15}. $$

I am fully aware that this is long overdue, but I wanted to suggest another way of finding the integral using a slightly different contour that has no indentations. This approach is quite similar to Brightsun's (basically the same!) but uses a little known fact about poles at the corner of a contour and how they should be treated. I have obviously skipped some intermediate steps in order to emphasize the different process. I hope someone will find this (at least partly) useful.

We'll start off by evaluating the following integral \begin{equation} \oint_C\frac{z^4}{-e^z-1}dz \end{equation} where $C$ is a rectangle of vertices $(0,-\pi),(R,-\pi),(R,\pi),(0,\pi)$.

Notice that the sum of the residues must be multiplied by $i\pi /2$ instead of $2\pi i$ since the two poles are at two corners of the contour of angle $\pi/2$ each. As a general rule, when the curve is smooth where it passes through a pole, the residue is multiplied by $i\pi$. However if the pole is at the tip of a corner of inner angle $\psi$ then the corresponding residue is multiplied by $i\psi $. The general form of the residue can be readily found:$$\text{Res}_{z=i(2k+1)\pi}\left(\frac{-z^4}{e^z+1}\right)=\left(\frac{-z^4}{e^z}\right)_{z=i(2k+1)\pi}=\left(\frac{-(i(2k+1)\pi)^4}{e^{i(2k+1)\pi}}\right)=(2k\pi+\pi)^4$$

The two poles on the left corners correspond to $k=0$ and $k=-1$ \begin{align} \int_0 ^R \frac{(x-i\pi)^4}{e^x - 1}\mathrm{d}x +& \int_{-\pi} ^{\pi} \frac{(R+iy)^4}{-e^{R+iy}- 1}i\mathrm{d}y +\\[8pt] &\int_{R}^0 \frac{(x+i\pi)^4}{e^{x}-1}\mathrm{d}x + \int_{\pi}^{-\pi}\frac{(iy)^4}{-e^{iy}-1}i\mathrm{d}y= i\frac{\pi}{2}\times\left(\pi^4+\pi^4\right) \end{align} \begin{align} \int_0 ^\infty \frac{(x-i\pi)^4-(x+i\pi)^4}{e^x - 1}\mathrm{d}x + \int_{\pi}^{-\pi}\frac{(iy)^4}{-e^{iy}-1}i\mathrm{d}y &=i \pi^5 \\[8pt]\end{align} As $R\to\infty$ the second integral vanishes. After manipulating the fourth integral to separate the imaginary from the real parts, we equate the imaginary parts of the two sides of the equation: \begin{align} \int_0 ^\infty \frac{(x-i\pi)^4-(x+i\pi)^4}{e^x - 1}\mathrm{d}x + \frac{1}{2} \int_{\pi}^{-\pi} y^4 \mathrm{d}y&= \pi^5 \\[8pt] -8\pi \int_0 ^\infty \frac{x^3-x\pi^2}{e^x - 1}\mathrm{d}x + \frac{1}{2} \int_{\pi}^{-\pi} y^4 \mathrm{d}y &= \pi^5\\[8pt] 8\pi I=\frac{4\pi^5}3+\frac {\pi^5}5-\pi^5=\pi^5\left(\frac{20+3-15}{15}\right)&=\frac{8\pi^5}{15} \end{align} $$\boxed{\int_0^\infty \frac{x^3}{e^x-1}dx=\frac{\pi^4}{15}}$$