[Python] Why can't sympy integrate this (Mathematica can)

Why can sympy not carry out this integral:

import sympy as sp
x = sp.symbols('x', real=True, nonzero=True)
sp.integrate(x**3/(sp.exp(x)-1), (x, 0, sp.oo))


this should return pi^4/15.

Continue reading...