An Integral involving $e^{ax} +1$ and $e^{bx} + 1$
$$I(a,b)=\int_{0}^{\infty}\frac{e^{ax}-e^{bx}}{x\left(e^{ax}+1\right)\left(e^{bx}+1\right)}dx = \int_{0}^{\infty} \left(\frac{1} {x\left(e^{bx}+1\right)}- \frac{1} {x\left(e^{ax}+1\right)} \right)dx$$
$$\frac{\partial I}{\partial a} = \int_0^{\infty} \frac{e^{ax}}{(1+e^{ax})^2} dx$$
Let $t = e^{ax}$. We get
$$\frac{\partial I}{\partial a} = \frac1{a} \int_1^{\infty} \frac{dt}{(1+t)^2} = \frac{1}{2a}$$
Similarly, $$\frac{\partial I}{\partial b} = - \frac1{b} \int_1^{\infty} \frac{dt}{(1+t)^2} = -\frac{1}{2b}$$
Hence, $I(a,b) = \frac1{2} \left( \log(a) - \log(b) \right) + c$.
Further $I(a,a) = 0$ and hence $c = 0$. Hence,
$$I(a,b) = \frac1{2} \log \left(\frac{a}{b}\right)$$
Here is a first solution, which uses no complex variable nor any differentiation, just plain real analysis. A second solution follows, which uses Fubini theorem.
For every positive $t$, let $I_t(a,b)$ denote the integral of the same function than in the integral which defines $I(a,b)$, but integrated from $t$ to $+\infty$. Then $$ I_t(a,b)=\int_{t}^{+\infty}\frac{\mathrm{d}x}{x\left(\mathrm{e}^{bx}+1\right)}-\int_{t}^{+\infty}\frac{\mathrm{d}x}{x\left(\mathrm{e}^{ax}+1\right)}= \int_{bt}^{at}\frac{\mathrm{d}x}{x\left(\mathrm{e}^{x}+1\right)}, $$ where the first decomposition of $I_t(a,b)$ into two parts is legal because we got rid of the problem at $x=0$ and the second formula stems from the changes of variables $x\leftarrow bx$ and $x\leftarrow ax$ in the two parts.
Extracting $1/(2x)$ from the last fraction and integrating, one gets $$ I_t(a,b)=\frac12\log\left(\frac{a}b\right)- \int_{bt}^{at}\frac{\mathrm{e}^{x}-1}{2x\left(\mathrm{e}^{x}+1\right)}\mathrm{d}x=\frac12\log\left(\frac{a}b\right)+O(t), $$ where the $O(t)$ term when $t\to0$ comes from the fact that the function in the last integral is bounded around $x=0$ and from the fact that this integral is over an interval of length $at-bt=O(t)$. Finally, when $t\to0$, $I_t(a,b)\to I(a,b)$ hence $$ I(a,b)=\frac12\log\left(\frac{a}b\right). $$
Here is a second solution, based on Fubini theorem. Introduce $$ u(c,x)=\frac{\mathrm{e}^{cx}}{(1+\mathrm{e}^{cx})^2}=\frac{\partial}{\partial c}\left(\frac{-1}{x(1+\mathrm{e}^{cx})}\right). $$ One sees that $$ \frac{\mathrm{e}^{ax}-\mathrm{e}^{bx}}{x\left(\mathrm{e}^{ax}+1\right)\left(\mathrm{e}^{bx}+1\right)}=\frac1{x\left(\mathrm{e}^{bx}+1\right)}-\frac1{x\left(\mathrm{e}^{ax}+1\right)}=\int_b^au(c,x)\mathrm{d}c, $$ hence $$ I(a,b)=\int_{0}^{\infty}\left(\int_b^au(c,x)\mathrm{d}c\right)\mathrm{d}x. $$ Since $u\ge0$, Fubini-Tonelli theorem says one can interchange the order of integration, hence $$ I(a,b)=\int_b^aU(c)\mathrm{d}c,\qquad U(c)=\int_{0}^{\infty}u(c,x)\mathrm{d}x. $$ Obviously, $$ u(c,x)=\frac{\partial}{\partial x}\left(\frac{-1}{c(1+\mathrm{e}^{cx})}\right),\qquad U(c)=\left.\frac{-1}{c(1+\mathrm{e}^{cx})}\right|^{x=+\infty}_{x=0}=\frac1{2c}, $$ and the value of $I(a,b)$ follows.
Here is my solution with the eta function:
Let $\text{Re}(s)>1$, and consider $$\int_{0}^{\infty}\frac{\left(e^{ax}-e^{bx}\right)x^{s-1}}{\left(e^{ax}+1\right)\left(e^{bx}+1\right)}dx=\int_{0}^{\infty}\frac{x^{s-1}}{e^{bx}+1}dx-\int_{0}^{\infty}\frac{x^{s-1}}{e^{ax}+1}.$$ Recall that for $\text{Re}(s)>1$ we have $$\Gamma(s)\eta(s)=\int_{0}^{\infty}\frac{t^{s-1}}{e^{t}+1}dt$$ where $\eta(s)$ is the dirichlet eta function. Hence, by a substitution, the above expression is $$\frac{\Gamma(s)\eta(s)}{b^{s}}-\frac{\Gamma(s)\eta(s)}{a^{s}}=\Gamma(s)\eta(s)\left(e^{-s\log b}-e^{-s\log a}\right)$$ $$=\Gamma(s)\eta(s)\left(s\left(\log a-\log b\right)+s^{2}\left(\cdots\right)+\cdots\right)$$ Since $\Gamma(s)$ has a simple pole with residue $1$ at $s=0$, and $\eta(0)=\left(1-2\right)\zeta(0)=\frac{1}{2}$, we have $$\lim_{s\rightarrow0}\int_{0}^{\infty}\frac{e^{ax}-e^{bx}}{\left(e^{ax}+1\right)\left(e^{bx}+1\right)}x^{s-1}dx=\frac{1}{2}\log\left(\frac{a}{b}\right).$$ Switching the limit and the integral is then justifiable using the dominated convergence theorem so that we conclude $$I(a,b)=\frac{1}{2}\log\left(\frac{a}{b}\right).$$