Integrating $\int_0^{\frac{\pi}{2}} x (\log\tan x)^{2n+1}\;dx$ Does anybody have any thoughts about how to integrate
$$I=\int_0^{\frac{\pi}{2}} x (\log\tan x)^{2n+1}\;dx$$
for integral $n$ where $n\ge 1$
In the case $n=0$
$$\int_0^{\frac{\pi}{2}} x \log\tan x \;dx=\lambda(3)=\frac{7}{8}\zeta(3)$$
I have managed to integrate the function when the exponent is even, that is $(\log\tan x)^{2n}$, using the substitution $y=\left(\frac{\pi}{2}-x \right)$ over the two intervals $[0,\frac{\pi}{4}]$ and $[\frac{\pi}{4},\frac{\pi}{2}]$, but the same trick does not apply in regard to the odd powers.
Basically via integration by parts I am left with the repeated integral
$$\int_0^{\frac{\pi}{2}} \int_0^x  (\log\tan u)^{2n+1}\;du\;dx$$
As far as I know $(\log\tan u)^{2n+1}$ does not have a definite integral I can use, so I am stuck. I've tried a few substitutions and those have not helped. Any ideas how I might proceed?
Some Added Background Notes


*

*To obtain a function more suitable for numerical integration use the substitution $u=\log \tan x$ to give


$$\int_0^{\frac{\pi}{2}} x (\log\tan x)^{n}\;dx=\int_{-\infty}^{+\infty} \arctan(e^u)\frac{u^n}{e^u+e^{-u}} \;du$$
This shows that the integral $I$ is closely related to the standard integral for the $\beta(n)$ function. The analogous integral $\int_0^{\infty} x (\log\tanh x)^{n}\;dx$ via a similar change of variables is seen to be related to the standard integral for the $\lambda(n)$ function.

By setting $x=\arctan u$ we are left with
$$ \mathcal{I}(n) = \int_{0}^{+\infty}\frac{\arctan u}{1+u^2}\left(\log u\right)^{2n+1}\,du =\left.\frac{d^{2n+1}}{d\alpha^{2n+1}}\int_{0}^{+\infty}\frac{u^\alpha \arctan u}{1+u^2}\,du\right|_{\alpha=0}$$
but the integral in the RHS is related to the Beta function. By un-doing the previous substitution,
$$ \int_{0}^{+\infty}\frac{u^\alpha \arctan u}{1+u^2}\,du = \int_{0}^{\pi/2}x\left(\sin x\right)^{\alpha}\left(\cos x\right)^{-\alpha}\,dx $$
where we may write $x$ as
$$ \arcsin(\sin x)= \sum_{n\geq 0}\frac{(\sin x)^{2n+1}}{(2n+1)4^n}\binom{2n}{n}$$
leading to:
$$ \int_{0}^{+\infty}\frac{u^{\alpha}\arctan u}{1+u^2}\,du=\sum_{n\geq 0}\frac{\binom{2n}{n}}{(2n+1)4^n}\cdot\frac{\Gamma\left(\frac{1}{2}-\frac{a}{2}\right)\, \Gamma\left(1+\frac{a}{2}+n\right)}{2\,\Gamma\left(\frac{3}{2}+n\right)}. $$
Now "it is enough" to differentiate both sides with respect to $\alpha$ the correct number of times and perform an evaluation at $\alpha=0$ in order to convert the original integral in a "twisted hypergeometric series", whose terms depend both on hypergeometric terms and generalized harmonic numbers (arising from the differentiation of the Beta function).
