I have solved the question. But I have problem with the condition on n.
This is what I think.
RHS converges if $n>1$(because of Zeta function).
For LHS to converge, the given function must tend to zero as $x$ tends to infinity. This is true for all $n\in R$. Also, the given function must converge to some value when $x$ tends to $0$. When $n=2$ the function on LHS tends to $1$ when $x$ tends to $0$. When $n>2$ it tends to zero. For $n<2$ it tends to infinity. So the basic conditions for convergence of integral on LHS are $n\geq 2$ (I have not checked convergence with any test since only this condition gives me the answer).
Why do the conditions on $n$ do not match on both sides of the equation?
Best Answer
In fact, the LHS integral is finite for $n=1+\epsilon$, $\forall\epsilon>0$. To see this, note that the integrand behaves like $\frac{1}{x^{1-\epsilon}}$ near $x=0$. Since $$ \int_0^1 \frac{1}{x^{1-\epsilon}}dx = \frac{1}{\epsilon}x^{\epsilon}\Big|^1_0=\frac{1}{\epsilon}<\infty, $$ by comparison test, we have $$\int_0^1 \frac{x^{1+\epsilon}e^x}{(e^x -1)^2}dx<\infty.$$ And since the integrand decays exponentially, we have $$ \int_0^\infty \frac{x^{1+\epsilon}e^x}{(e^x -1)^2}dx<\infty,\quad\forall\epsilon>0. $$ This shows that the both sides converge for $n>1$.