I have been trying to find the value of:$$\int_0^\infty\frac{\tanh(ax)}{e^x-1}dx=\int_0^\infty\frac{e^{2ax}-1}{(e^{2ax}+1)(e^x-1)}dx$$ Under u-substitution:
Let $u=e^{ax}$, $x=\frac{\log(u)}{a}$, $dx=\frac{du}{ua}$
$$\int_1^\infty \frac{u^2-1}{(u^2+1)(u^{\frac{1}{a}}-1)}\frac{du}{ua}$$
I don't really know where to proceed further from here or if I should use a series expansion of $\tanh(ax)$.
I would rewrite it as such:
$$\int_0^\infty \tanh (a x)\frac{e^{-x}}{1-e^{-x}}dx=\sum_{n=1}^\infty\int_0^\infty \tanh ax \,e^{-nx}dx$$ $$=\sum_{n=1}^\infty\frac1a\int_0^\infty \tanh u \,e^{-nu/a}du=\sum_{n=1}^\infty \frac1 a f\left(\frac{n}{a}\right)$$ Now we need to solve for a different function where Mathematica helps us to find the solution in terms of special digamma function $\psi(x)=\frac{\Gamma'(x)}{\Gamma(x)}$: $$f(x)=\int_0^\infty \tanh u\, e^{-xu}du=-\frac1x+\frac12\left[\psi\left(\frac{x}{4}+\frac12\right)-\psi\left(\frac{x}{4}\right)\right]$$ It looks bad, but fortunately, we have $\psi(x+1)=\psi(x)+\frac1x$ which will simplify things in a sum for nice values of $a$.
The overall sum: $$\sum_{n=1}^\infty \frac1 a f\left(\frac{n}{a}\right)=\sum_{n=1}^\infty\left[-\frac{1}{n}+\frac1{2a}\left[\psi\left(\frac{n}{4a}+\frac12\right)-\psi\left(\frac{n}{4a}\right)\right]\right]$$
If $a$ is half integer: $a=q/2$, rewrite $n=qm+p$ (summing now over all $m\geq 0$ and $0<p\leq q$): $$\sum_{n=1}^\infty\left[-\frac{1}{n}+\frac1{q}\left[\psi\left(\frac{qm+p}{2q}+\frac12\right)-\psi\left(\frac{qm+p}{2q}\right)\right]\right]=$$ $$=\frac1{q}\sum_{n=1}^\infty\left[-\frac{q}{n}+\left[\psi\left(\frac{m+1}{2}+\frac{p}{2q}\right)-\psi\left(\frac{m}{2}+\frac{p}{2q}\right)\right]\right]=$$ $$=\frac1{q}\sum_{n=1}^\infty\left[-\frac{q}{n}+\left[\psi\left(\frac{m+1}{2}+\frac{p}{2q}\right)-\psi\left(\frac{m}{2}+\frac{p}{2q}\right)\right]\right]=$$ $$=\frac1{q}\sum_{m=0}^\infty\sum_{p=1}^q\left[-\frac{q}{qm+p}+\left[\psi\left(\frac{m+1}{2}+\frac{p}{2q}\right)-\psi\left(\frac{m}{2}+\frac{p}{2q}\right)\right]\right]$$
We cannot split off the first part of the sum because the harmonic sum diverges. But we can exchange the order, and then group two adjacent terms, $m=2k+0$ and $m=2k+1$: $$=\frac1{q}\sum_{p=1}^q \sum_{m=0}^\infty\left[-\frac{q}{qm+p}+\left[\psi\left(\frac{m+1}{2}+\frac{p}{2q}\right)-\psi\left(\frac{m}{2}+\frac{p}{2q}\right)\right]\right]= \quad(*)$$ $$=\frac1{q}\sum_{p=1}^q \sum_{k=0}^\infty\left[-\frac{q}{2kq+p}-\frac{q}{(2k+1)q+p}+\left[\psi\left(k+1+\frac{p}{2q}\right)-\psi\left(k+\frac{p}{2q}\right)\right]\right]=$$ I used the telescoping property on the $\psi$ terms. Now we recognize the inner bracket as $1/(k+p/(2q))$ from the properties of the digamma function. $$=\frac1{q}\sum_{p=1}^q \sum_{k=0}^\infty\left[-\frac{q}{2kq+p}-\frac{q}{(2k+1)q+p}+\frac{2q}{2kq+p}\right]=$$ $$=\sum_{p=1}^q \sum_{k=0}^\infty\left[\frac{1}{2kq+p}-\frac{1}{(2k+1)q+p}\right]=\sum_{p=1}^q \sum_{m=0}^\infty\frac{(-1)^m}{mq+p}=\sum_{n=1}^\infty \frac{b_n}{n}$$ where $b_n$ is a sequence of $\pm 1$: $$(b_n)=\{\underbrace{1,1,1,1,1,\ldots}_q,\underbrace{-1,-1,-1,-1,-1,\ldots}_q,\ldots\}$$ This is a Dirichlet series of sequence $b_n$, evaluated at $1$ and can be evaluated for special cases of $q$ quite efficiently (such as the $\ln 2$ from alternating harmonic series at $q=1$ and all the solutions stated by @Claude), but I think this form is especially elegant.
$$\int_0^\infty \frac{\tanh a x}{e^x-1}dx=\sum_{n=1}^\infty\frac{b^{(2a)}_n}{n};\quad 2a\in \mathbb{N}$$
This derivation was probably "around the corner" and much more complex than necessary. There are easier ways, one of which can be seen from the power series of the integrand under $u=e^{-x}$ substitution.
A "closed form" can be computed for the integer case if we resume the calculation $(*)$ and go in a different direction, writing the infinite sum over $m$ as a limit of a finite sum, using $\sum_{m=0}^M \frac{1}{m+p/q}=\psi(M+1+p/q)-\psi(p/q)$ and recognizing the telescoping sum in the last bracket of $(*)$:
$$\frac1{q}\sum_{p=1}^q \lim_{M\to\infty}\left[-\psi\left(M+1+\frac{p}{q}\right)+\psi\left(\frac{p}{q}\right)+\psi\left(\frac{M+1}{2}+\frac{p}{2q}\right)-\psi\left(\frac{p}{2q}\right)\right]$$
Asymptotic expansion $\psi(x)\approx \ln x + \mathcal{O}(1/x)$ helps us derive the exact value of the limit above:
$$\frac1{q}\sum_{p=1}^q \left[-\ln 2+\psi\left(\frac{p}{q}\right)-\psi\left(\frac{p}{2q}\right)\right]$$
Now we can take out $\ln 2$, use known finite sums for rational arguments of $\psi$:
$$-(\gamma+\ln q)-\ln 2-\frac1{q}\sum_{p=1}^q \left[\psi\left(\frac{p}{2q}\right)\right]$$ and the leftover term can be processed by taking advantage of the reflection formula $\psi(1-x)-\psi(x)=\pi \cot\pi x$. Consider the sum (same formula as before):
$$\sum_{p=1}^{2q} \psi\left(\frac{p}{2q}\right)=-2 q (\gamma + \ln 2q)$$ $$\sum_{p=1}^{q} \left[\psi\left(\frac{p}{2q}\right)+\psi\left(1-\frac{p}{2q}\right)\right]+\psi(1)-\psi(1/2)=-2 q (\gamma + \ln 2q)$$ $$\sum_{p=1}^{q} \left[2\psi\left(\frac{p}{2q}\right)+\pi \cot \pi \frac{p}{2q}\right]-\gamma+2\ln 2+\gamma=-2 q (\gamma + \ln 2q)$$ $$\sum_{p=1}^{q} \psi\left(\frac{p}{2q}\right)=- q (\gamma + \ln 2q)-\ln 2-\frac{\pi}{2}\sum_{p=1}^q \cot \frac{\pi p}{2q}$$ The final closed form of the original sum for half-integer arguments is therefore simply: $$-(\gamma+\ln q)-\ln 2+\frac1{q}\left[q (\gamma + \ln 2q)+\ln 2+\frac{\pi}{2}\sum_{p=1}^q \cot \frac{\pi p}{2q}\right]=$$ $$=\frac1{q}\left[\ln 2+\frac{\pi}{2}\sum_{p=1}^q \cot \frac{\pi p}{2q}\right]$$
$$\int_0^\infty \frac{\tanh a x}{e^x-1}dx=\frac1{2a}\left[\ln 2+\frac{\pi}{2}\sum_{p=1}^{2a} \cot \frac{\pi p}{4a}\right]$$
This reproduces the table of @Claude in a closed form. The sum runs across cotangents of evenly divided $\pi/2$ interval into $2a$ segments, and is a Riemann sum of $\int_0^{\pi/2}\cot x\,dx$ which diverges. This term is the one that causes divergent behaviour for large $a$ (logarithmic divergence).
For arbitrary real-valued $a$, I'm almost certain that there is no simple closed form. But an asymptotic expansion can be obtained through Euler-Maclaurin's formula, truncated to a reasonable number of terms, for example,
$$\int_0^\infty \frac{\tanh a x}{e^x-1}dx\asymp \frac{1}{2a}\ln 2 -\ln\sin\frac{\pi}{4a}+\frac12\frac{\pi}{4a}\cot\frac{\pi}{4a} +\frac1{12}\left(\frac{\pi}{4a}\cot\frac{\pi}{4a}\right)^2 $$
which is accurate to about a percent.
Not an answer.
Since marty cohen started with some values, I shall add a few other (including those already given) $$\left( \begin{array}{cc} a & \text{result} \\ \frac{1}{2} & \log (2) \\ 1 & \frac{1}{4} (\pi +2\log (2)) \\ \frac{3}{2} & \frac{1}{9} \left(2 \sqrt{3} \pi +3\log (2)\right) \\ 2 & \frac{1}{8} \left(\pi +2 \sqrt{2} \pi +2\log (2)\right) \\ 3 & \frac{1}{36} \left(\left(15+4 \sqrt{3}\right) \pi +6\log (2)\right) \\ 4 & \frac{1}{16} \left(\pi +2 \sqrt{2} \pi +8 \pi \cos \left(\frac{\pi }{8}\right)+2\log (2)\right) \end{array} \right)$$ Using numerical integration (let $a=10^k$), we have as results $$\left( \begin{array}{cc} k & \text{integral} \\ 0 & 1.13197 \\ 1 & 3.15568 \\ 2 & 5.42741 \\ 3 & 7.72688 \\ 4 & 10.0292 \\ 5 & 12.3317 \\ 6 & 14.6343 \\ 7 & 16.9369 \\ 8 & 19.2395 \\ 9 & 21.5420 \\ 10 & 23.8446 \end{array} \right)$$ which seems to be almost a straight line when plotted as a function of $k$ $(R^2=0.99997)$ and the slope seems to be very close to $\log(10)$. $$\begin{array}{clclclclc} \text{} & \text{Estimate} & \text{Standard Error} & \text{Confidence Interval} \\ a & 0.928645 & 0.049392 & \{0.814746,1.042544\} \\ b & 2.287000 & 0.008349 & \{2.267748,2.306253\} \\ \end{array}$$
