It is well known that for $n\in\mathbb{N}$ and $n>0$ (an maybe even for more than these restrictions): $$I_n = \int_0^\infty\frac{x^n}{e^x-1}dx = \zeta(n+1)n!$$ which, analytically can be shown easily by expanding $1/(1-e^{-x})$ into a geometric series, which leads to trivial integrals, and by using $\zeta(n+1)=\sum_{l=1}^\infty l^{-(n+1)}$. So, eg.: $$I_1 = \pi^2/6$$ $$I_2 = 2\zeta(3)$$ a.s.o...
Now, if I try even the simplest case with sage, I get this 'nifty' little results
sage: integrate(x/(exp(x)-1),x,0,oo)
-1/6*pi^2 + limit(-1/2*x^2 + x*log(-e^x + 1) + polylog(2, e^x), x,
+Infinity, minus)
Is there any trick to simplify this down to the final result, or is this about as far as I can get with sage alone?
PS.: it is probably needless to say that (once again ... :( ...)
In[1]:= Integrate[x/(Exp[x] - 1), {x, 0, Infinity}]
Out[1]:= Pi^2/6
