Here is a more elementary solution:
Rewrite the given integral as:
$$I=\int_0^{\infty} \frac{x^2e^{-x}}{1+e^{-x}}\,dx$$
Use the substitution ##e^{-x}=t## to obtain:
$$I=\int_0^1 \frac{\ln^2 t}{1+t}\,dt$$
Next, use integration by parts in the following way:
$$I=\left(\ln^2t \ln(1+t)\right|_0^1-2\int_0^1 \frac{2\ln t\ln(1+t)}{t}\,dt$$
Notice that the first term is zero, so we are left with:
$$I=-2\int_0^1 \frac{\ln t \ln(1+t)}{t}\,dt$$
Since,
$$\ln(1+t)=-\sum_{k=1}^{\infty} \frac{(-1)^kt^k}{k}$$
Hence,
$$I=2\sum_{k=1}^{\infty} \frac{(-1)^k}{k}\int_0^1 \frac{t^k\ln t}{t}\,dt$$
It is not difficult to show that
$$\int_0^1 \frac{t^k\ln t}{t}\,dt =-\frac{1}{k^2}$$
From above, we get:
$$I=-2\sum_{k=1}^{\infty} \frac{(-1)^k}{k^3}$$
It can be easily shown that
$$\sum_{k=1}^{\infty} \frac{(-1)^k}{k^3}=-\frac{3}{4}\zeta(3)$$
Hence, the final answer is:
$$I=\frac{3}{2}\zeta(3)$$
which agrees with post #2.
I hope that helps.