A closed form for an integral expressed as a finite series of $\zeta(2k+1)$, $\pi^m$ and a rational?

I think the following articles can give a clue: http://www.tandfonline.com/doi/abs/10.1080/10652460701688125?journalCode=gitr20 (Closed-form evaluation of some families of cotangent and cosecant integrals, by D. Cvijović) and http://www.sciencedirect.com/science/article/pii/S0377042702003588 (Integral representations of the Riemann zeta function for odd-integer arguments, by D. Cvijović and J. Klinowski).

Update. In particular, using $2x^n=E_n(x)+E_n(1+x)$ and $E_n(1+x)=\sum\limits_{k=0}^n\binom{n}{k}E_k(x)$, we get $$f(n)=\frac{\pi^{n+2}}{2}\left[I_n+\sum\limits_{k=0}^n\binom{n}{k}I_k-I_{n+1}-\sum\limits_{k=0}^{n+1}\binom{n+1}{k}I_k\right]=\frac{\pi^{n+2}}{2}\left[I_n-2I_{n+1}-\sum\limits_{k=1}^n\binom{n}{k-1}I_k\right], \tag{1}$$ where $$I_n=\int\limits_0^1E_n(x)\csc{(\pi x)}\,dx.$$ Odd indices in (1) doesn't contribute because $E_n(1-x)=(-1)^nE_n(x)$ leads to the integrand which is antisymmetric with regard to $x\to 1-x$. For even indices, $I_n$ was calculated in the second paper indicated above as $$I_{2n}=(-1)^n\frac{(4-2^{1-2n})(2n)!}{\pi^{2n+1}}\zeta(2n+1),$$ and we obtain, for even $n=2m$: $$f(2m)=(-1)^m\frac{\pi}{2}(4-2^{1-2m})(2m)!\,\zeta(2m+1)- \sum\limits_{p=1}^m(-1)^p\;\frac{\pi^{2(m-p)+1}}{2}\;\binom{2m}{2p-1}(4-2^{1-2p})(2p)!\,\zeta(2p+1),$$ and for odd $n=2m-1, m>1$: $$f(2m-1)=(-1)^{m+1}(4-2^{1-2m})(2m)!\,\zeta(2m+1)-\sum\limits_{p=1}^{m-1}(-1)^p\;\frac{\pi^{2(m-p)}}{2}\;\binom{2m-1}{2p-1}(4-2^{1-2p})(2p)!\,\zeta(2p+1).$$

Not an answer, but... Mathematica evaluates the indefinite integral $$\int x^k \csc x d x$$ in terms of polylogs, from which your integral follows readily, so the obvious question is - how does it do the indefinite integral. The short answer is: The Risch algorithm (or, more precisely, an extension, since the original algorithm only deals with elementary antiderivatives). However, if you want to do it by hand, the two methods which come to mind are:

1. Make the substitution $u = \exp(i x),$ which transforms your integrand into a power of log times a rational function - this is plausbible, since the arguments of the polylogs have the form $\exp(i x),$ so Mathematica might, indeed, be doing this.
2. Expand $\csc x$ in its Laurent series around $0,$ and integrate term by term. The Laurent series is explicit (in terms of Bernoulli numbers), so it may be that a simple manipulation gets it to your linear combination of zetas (I haven't tried).