Denote your integral as $\mathfrak{I}(n)$ and apply IBP by choosing $u=x^{2n+1}$ and $\mathrm dv=\cot(x)\mathrm dx$ to get
\begin{align*}
\mathfrak{I}(n)&=\int_0^{\pi/2}x^{2n+1}\cot(x)\mathrm dx =\underbrace{\left[(2n+1)\cdot x^{2n}\log(\sin x)\right]_0^{\pi/2}}_{\to0}-(2n+1)\int_0^{\pi/2}x^{2n}\log(\sin x)\mathrm dx\\
&=-(2n+1)\int_0^{\pi/2}x^{2n}\log(\sin x)\mathrm dx
\end{align*}
Now utilizing the well-known Fourier series expansion of $\log(\sin x)$, which converges within $[0,\pi]$, and switching the order of summation and integration further gives us
\begin{align*}
\mathfrak{I}(n)&=-(2n+1)\int_0^{\pi/2}x^{2n}\log(\sin x)\mathrm dx\\
&=-(2n+1)\int_0^{\pi/2}x^{2n}\left[-\log(2)-\sum_{k=1}^\infty\frac{\cos(2kx)}k\right]\mathrm dx\\
&=\log(2)\left(\frac\pi2\right)^{2n+1}+(2n+1)\sum_{k=1}^\infty\frac1k\underbrace{\int_0^{\pi/2}x^{2n}\cos(2kx)\mathrm dx}_{=J}
\end{align*}
The integral $J$ can be computed via IBP again which explains the connection to values of the Riemann Zeta Function hence for integer $n$ every IBP step produces another reciprocal power of $n$ which overall combines to sums that can be expressed with the help of the Riemann Zeta Function.
As one may see the values for $n=0$ and $n=1$ can be easily verfied since for $n=0$ $J$ is overall $0$ aswell whereas for $n=1$ the latter integral can be expressed using the Dirichlet Eta Function. To be precise we got
\begin{align*}
n=0:~~~\mathfrak{I}(0)&=\log(2)\left(\frac\pi2\right)^{1}+(1)\sum_{k=1}^\infty\frac1k\underbrace{\int_0^{\pi/2}\cos(2kx)\mathrm dx}_{=0}\\
&=\frac{\pi\log(2)}2
\end{align*}
\begin{align*}
n=1:~~~\mathfrak{I}(1)&=\log(2)\left(\frac\pi2\right)^{3}+(2+1)\sum_{k=1}^\infty\frac1k\int_0^{\pi/2}x^2\cos(2kx)\mathrm dx\\
&=\log(2)\left(\frac\pi2\right)^{3}+3\sum_{k=1}^\infty\frac1k\left[\frac\pi4\frac{\cos(\pi k)}{k^2}\right]_0^{\pi/2}\\
&=\log(2)\left(\frac\pi2\right)^{3}-\frac{3\pi}4\sum_{k=1}^\infty\frac{(-1)^{k+1}}{k^3}\\
&=\log(2)\left(\frac\pi2\right)^{3}-\frac{3\pi}4\eta(3)\\
&=\log(2)\left(\frac\pi2\right)^{3}-\frac{9\pi}{16}\zeta(3)\\
&=\frac1{16}(\pi^3\log(4)-9\pi\zeta(3))
\end{align*}
Note that we used the relation $\eta(s)=(1-2^{1-s})\zeta(s)$. Similiar can be done for all integer $n$. So as at least close to a closed-form I can offer the following formula
$$\therefore~\mathfrak{I}(n)~=~\log(2)\left(\frac\pi2\right)^{2n+1}+(2n+1)\sum_{k=1}^\infty\frac1k\int_0^{\pi/2}x^{2n}\cos(2kx)\mathrm dx$$