How to find $f(c)=\int_{0}^{\infty}\frac{x^2}{e^{cx}-e^x}dx$? I am trying to find out how $$ f(c)=\int_{0}^{\infty}\frac{x^2}{e^{cx}-e^x}dx $$ depends on $c$, for $c>1$. I haven't been able to compute the antiderivative of the integrand. Does anybody see how to compute the improper integral? 
Thank you! 
 A: $\newcommand{\bbx}[1]{\bbox[8px,border:1px groove navy]{{#1}}\ }
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{\displaystyle{#1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\mc}[1]{\mathcal{#1}}
 \newcommand{\mrm}[1]{\mathrm{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
\left.\int_{0}^{\infty}{x^{2} \over \expo{cx} - \expo{x}}\,\dd x
\,\right\vert_{\ c\ >\ 1} & =
\int_{0}^{\infty}{x^{2}\expo{-cx} \over 1 - \expo{-\pars{c - 1}x}}\,\dd x
\end{align}
With the substitution
$\ds{t \equiv \expo{-\pars{c - 1}x} \iff x = -\,{\ln\pars{t} \over c - 1}}$:
\begin{align}
\left.\int_{0}^{\infty}{x^{2} \over \expo{cx} - \expo{x}}\,\dd x
\,\right\vert_{\ c\ >\ 1} & =
\int_{1}^{0}{\bracks{-\ln\pars{t}/\pars{c - 1}}^{\, 2}\, t^{c/\pars{c - 1}}
\over 1 - t}\,\bracks{-\,{\dd t \over \pars{c - 1}t}}
\\[5mm] & =
{1 \over \pars{c - 1}^{3}}
\int_{0}^{1}{t^{1/\pars{c - 1}}\ln^{2}\pars{t} \over 1 - t}\,\dd t
\\[5mm] & =
\left. -\,{1 \over \pars{c - 1}^{3}}\,\partiald[2]{}{\mu}
\int_{0}^{1}{1 - t^{\mu} \over 1 - t}\,\dd t
\,\right\vert_{\ \mu\ =\ 1/\pars{c - 1}}\label{1}\tag{1}
\\[5mm] & =
-\,{1 \over \pars{c - 1}^{3}}\,\Psi\, ''\pars{{1 \over c - 1} + 1}
\\[5mm] & =\
\bbx{-\,{1 \over \pars{c - 1}^{3}}\,\Psi\, ''\pars{c \over c - 1}}
\end{align}

$\ds{\Psi\ \mbox{is the}\ Digamma\ Function}$. Note the well known identity
  $\ds{\left.\Psi\pars{\mu + 1} + \gamma =
\int_{0}^{1}{1 - t^{\mu} \over 1 - t}\,\dd t
\,\right\vert_{\ \Re\pars{\mu}\ >\ -1}}$ in expression \eqref{1}. $\ds{\gamma}$ is the Euler-Mascheroni Constant.

A: Since $c>1$, we may write the integrand function as a geometric series:
$$ \frac{x^2}{e^{cx}-e^{x}}=x^2 e^{-x}\sum_{n\geq 1}e^{n(1-c)x} \tag{1}$$
then perform termise integration:
$$ \int_{0}^{+\infty}\frac{x^2}{e^{cx}-e^x}\,dx = \color{red}{\sum_{n\geq 1}\frac{2}{(1+(c-1)n)^3}}.\tag{2} $$
In particular, for $c=2$ we get $2\,\zeta(3)-2$, for $c=3$ we get $\frac{7}{4}\zeta(3)-2$ and for $c=4$ we get $\frac{26}{27}\zeta(3)+\frac{4\pi^3}{81\sqrt{3}}-2$. Additionally, for integer values of $c$ our function $f$ is related with the trilogarithm (inheriting peculiar reflection and duplication formulas) through a discrete Fourier transform.
