Well! This integral is really interesting and has a lot to do with Bernoulli's Numbers, Riemann Zeta function, Polygamma function, in the last what you were trying is nothing less than Laplace Transformation
$$I = \int_0^{\infty}\frac {t}{e^t-1}dt$$
First I'll try to answer why $$I =-\frac{\partial}{\partial s}\int_0^{\infty}\frac {e^{-st}}{e^t-1}dt$$ is failing to conclude answer and after that, I'll try to add few other ways you can conclude to the answer.
- Why $I =-\frac{\partial}{\partial s}\int_0^{\infty}\frac {e^{-st}}{e^t-1}dt$ fails
Answer On the very first look you can conclude that the definite integral $\int_0^{\infty}\frac {e^{-st}}{e^t-1}dt$ doesn't exist(Undefined)
Hint Lower limit of definite integral is $0$ and when $t =0 $ numerator = $1$ but denominator $=0$ However, for $\int_0^{\infty}\frac {t}{e^t-1}dt$ it's totally different
Even if you try to find by taking laplace transformation you'll eventually end up for $\Gamma(0)$
How to solve using Transformation only?
This may answer your question but as you asked for a few other ways to solve this integral I'm adding a few transformation based answers.
$1$] Using Laplace Transform Identity
I'm unable to provide a link from where I found this but if I could remember you can find it in "Handbook of Mathematics for Engineer's & Scientists" by Andrei & Alexander
If
$$\begin{align*}
F(s) = \mathcal{L}f(t)
&= \int_0^{\infty}e^{-st}f(t)dt\\
& \implies \sum_{s=1}^{\infty}F(s) = \int_0^{\infty} \frac {f(t)}{e^t-1}dt\\
& \text{We can put } f(t) = t \implies \mathcal{L}(t) = \frac 1{s^2}\\
& \text{Thus,}\\
& \color{green}{I = \int_0^{\infty} \frac{t}{e^t-1}dt = \sum_{s=1}^{\infty}\frac {1}{s^2}=\zeta(2) = \frac {\pi^2}{6}}\\
\end{align*}$$
$$\begin{align*}
I
&= \int_0^{\infty}\frac t{e^t-1}dt\\
& \text{We substitute, } e^{-t} = x \implies t = -\ln x \implies dt = -\frac{dx}{x} \\
& \color{blue}{I = -\int_0^1 \frac {\ln x}{1-x}dx}\\
\end{align*}$$
Now,
$$\color{red}{\psi_0 (z) = -\gamma + \int_0^{1} \frac {1-x^{z-1}}{1-x}dx}$$
We'll differentiate and that is trigamma function $\psi_1(z)$
$$\implies \frac {\partial\psi_0}{\partial z}= \psi_1 (z) = -\int_0^1 \frac {x^{z-1}\ln x}{1-x}dx$$
$$\color{green}{I = \left[\psi_1(z)\right]_{z=1} = \frac {\pi^2}{6}}$$
- $3]$ Using the integral form of even values of zeta function
$$\left|\zeta(2n) = \frac 1{(2n-1)!}\int_0^{\infty}\frac{x^{2n-1}}{e^x-1}dx\right|_{n =1}$$
- $4]$ More general form of the above integral in terms of Bernoulli's Number
$$\left|\int_0^{\infty} \frac {x^{2n-1}}{e^{px}-1}dx = (-1)^{n-1}\left(\frac {2\pi}{p}\right)^{2n} \frac{B_{2n}}{4n}\right|_{n, p = 1, 1}$$
$B_2 = \frac 1{6}$
- My last attempt: The cosine Transformation
$$f_c(u) = \int_0^{\infty}f(x)\cos(ux)dx$$
for $$\text{if }f(x) = \frac x{e^{ax}-1} \implies f_c(u) = \frac 1{2u^2} - \frac {\pi^2}{2a^2\sinh^2(\pi a^{-1}u)} $$
I tried to take the limit as $u \to 0$ but didn't work maybe because cosine transformation is defined for $u \in (0, \infty)$ but at $u = 0.0001$ it's $f_c(u = 0.0001) = 1.644934043288231 ≈ \pi^2/6$
I'm outsourcing most of the work from another answer. Perform integration by parts with $f = \ln(\frac{1}{x}-1) =\ln(1-x)-\ln(x) $ and $g' = \frac{\ln(1+x)}{x}$ to observe
$$\begin{align}\int_0^{1/2} \frac{1}{x}\ln(1+x)\ln\left(\frac{1}{x}-1\right)dx &= -\ln\left(\frac{1}{x}-1\right)\text{Li}_2(-x)\bigg\vert_{0}^{1/2} -\int_0^{1/2}\frac{\text{Li}_2(-x)}{x}+\frac{\text{Li}_2(-x)}{1-x} dx\\
&=-\int_0^{1/2}\frac{\text{Li}_2(-x)}{x}dx-\int_0^{1/2}\frac{\text{Li}_2(-x)}{1-x} dx\\ \end{align}$$
We have $$-\int_0^{1/2}\frac{\text{Li}_2(-x)}{x}dx =-\text{Li}_3\left(-\frac{1}{2}\right)$$ more or less by definition and the second is seen here to be given by
$$-\int_0^{1/2}\frac{\text{Li}_2(-x)}{1-x}dx = \text{Li}_3\left(-\frac{1}{2}\right) + \frac{13}{24}\zeta(3) $$
Edit: Looking at Bob's answer, it is really the calculation of his integral that is most of the work in the separate question I linked.
Best Answer
Indeed, the issue you're running into with your approach is that you have taken a convergent integral and split it up into two divergent integrals.
For another approach we first define $$ I=\int_0^1\frac{x^a-1}{\log x}\,\mathrm dx. $$ Expanding $x^a=e^{a\log x}$ yields $$ I=\sum_{k=1}^\infty\frac{a^k}{k!}\int_0^1\log^{k-1}x\,\mathrm dx=-\sum_{k=1}^\infty\frac{(-a)^k}{k!}\int_0^1(-\log x)^{k-1}\,\mathrm dx. $$ A substitution of $u=-\log x$ turns the remaining integral into a gamma function giving $$ I=-\sum_{k=1}^\infty\frac{(-a)^k}{k!}(k-1)!=-\sum_{k=1}^\infty\frac{(-a)^k}{k}=\log(1+a), $$ which follows form the series expansion of the logarithm.