I would suggest splitting up x³/(e^x-1) as:

x³ e^-x/(1-e^-x)=x³(e^-x+e^-2x+e^-3x+...)

i.e. a geometric series (which we can do since e^-x < 1 for x > 0).

Note that ò₀^¥ x³ e^{-n x} dx=ò₀^¥ u³ e^-u/n⁴ du = 6/n⁴.

I think it is easy if you're willing to accept Euler's sine product. If we start with:

sin x=x(1-x²/p²)(1-x²/(4p²))...

then on taking logarithms and differentiating you get:

cot x=1/x+

å

[1/(x+rp)+1/(x-rp)]

where r is taken over the naturals.

Differentiate each side 2n-1 times:

d2n-1/dx2n-1 cot x=-(2n-1)!(1/x2n+

å

1/(x+rp)2n+

å

1/(x-rp)2n)

Bring the 1/x²ⁿ term to the other side then take the limit as x® 0:

lim x® 0

((d2n-1/dx2n-1 cot x)+(2n-1)!/x2n)=-2(2n-1)!z(2n)/ p2n

Write cot x=i(e^{2i x}+1)/(e^{2i x}-1), multiply throughout by -i and you get:

lim x® 0

(d2n-1/dx2n-1((e2i x+1)/(e2i x-1))-i(2n-1)!/x2n) = 2i(2n-1)!z(2n)

So:

lim x® 0

(d2n-1/dx2n-1((e2i x+1)/(2i x)×2i x/(e2i x-1)) -i(2n-1)!/x2n)=2i(2n-1)!z(2n)

Setting y=2i x and using the chain rule this becomes:

lim y® 0

(d2n-1/dy2n-1((ey +1)/y×y/(ey-1))×22n-1 (-1)n+1i+22n(-1)n+1(2n-1)!i/y2n)=2i(2n-1)!z(2n)

(*)

Note that (e^y+1)/y and y/(e^y-1) can both be series expanded:

(e^y+1)/y=2/y+1+y/2!+y²/3!+y³/4!

y/(e^y-1)=B₀+B₁ y+B₂ y²/2!+...

where B_n is the nth Bernoulli number (using the well known fact that y/(e^y-1) is the generating function for the Bernoulli numbers).

If you substitute these series in (*) and expand out the product then you should find that the reciprocal terms cancel (so luckily the function doesn't blow up at the origin) and then the final result follows on applying the (well known??) identity:

n! Bn=

n å r=1

C(n,r)Br

which is easily proved by equating coefficients in the identity:

x/(e^x-1)×e^x=x+x/(e^x-1)

(e^x+1)/(e^2x-1)=1/(e^x-1).