In some physics problem I ran into the integral with a non-negative integer $n$ of the form: \begin{align} I(n) \equiv \int_0^\infty d x\, x^{2n+1} \tanh(\pi x)\, \log (1 - e^{-2\pi x}) \ , \end{align} which I cannot perform even with Mathematica for general $n$, but for small $n$ the values are given analytically (e.g. in Eq. (3.10) and (3.23) of https://arxiv.org/pdf/1708.00305.pdf): \begin{align} \begin{aligned} I(0) &= - \frac{1}{8}\,\log 2 + \frac{9}{16\pi^2}\,\zeta(3) \ ,\\ I(1) &= - \frac{1}{64}\,\log 2 - \frac{3}{32\pi^2}\,\zeta(3) + \frac{225}{128\pi^4}\,\zeta(5) \ , \end{aligned} \end{align}
It is likely that the integral generally takes the form \begin{align} I(n) = c_0\,\log 2 + \sum_{i=1}^{n+1}\,c_i\,\zeta(2i+1)\ . \end{align} Can we fix the values of the coefficients $c_i~(i=0,1,\cdots,n+1)$ for integer $n\ge 0$?
Thank you for your help in advance.