Formula for the general case of $\sum_{n=1}^\infty \frac{H_n}{n^2} = 2\zeta(3)$

Question

Introduction

I'm reading "Advanced Integration Techniques" by Zaid Alyafeai, and he brings up quite a interesting example regarding Euler sums. The definition is different in the book and on WolramMathWorld, but in the book he defines Euler sums as infinite sums of the form:

$$S_p(r,n) = \sum_{k=1}^\infty \frac{(H_k^{(p)})^r}{k^n}$$

where $H_k^{(p)}$ is a generalized harmonic number. This is quite a general case so let's look at some more specific cases. Let $\vartheta(n)$

$$\vartheta(n) = S_1(1,n) = \sum_{k=1}^\infty \frac{(H_k^{(1)})^1}{k^q} = \sum_{k=1}^\infty \frac{H_k}{k^n}$$

The interesting example I talked about was that

$$\vartheta(2) = \sum_{k=1}^\infty \frac{H_k}{k^2} = 2\zeta(3)$$

I'm going to quickly derive this so you get the idea. Start by using the integral representation of the Harmonic numbers, so we get:

\begin{align} \sum_{k=1}^\infty \frac{H_k}{k^2} & = \int_0^1 \frac{1}{1-x}\sum_{k=1}^\infty \frac{1-x^k}{k^2} dx \\\\ & = \int_0^1 \frac{\zeta(2)-\operatorname{Li}_2(x)}{1-x} dx\end{align}

Here $\operatorname{Li}_2(x)$ is the Polylogarithm of second order. We now use one of the functional equations for $\operatorname{Li}_2(x)$ which states that

$$\zeta(2)-\operatorname{Li}_2(x) = \operatorname{Li}_2(1-x)+\ln(x)\ln(1-x)$$

Plugging that into the integral, we get

\begin{align} \int_0^1 \frac{\zeta(2)-\operatorname{Li}_2(x)}{1-x} dx & = \int_0^1 \frac{ \operatorname{Li}_2(1-x)+\ln(x)\ln(1-x)}{1-x} \\\\ & = \int_0^1 \frac{\operatorname{Li}_2(1-x)}{1-x} dx + \int_0^1 \frac{\ln(x)\ln(1-x)}{1-x} dx \\\\ & = \zeta(3) + \zeta(3) = 2\zeta(3)\end{align}

Now, because of this beautiful result, I started wondering what e.g. $\vartheta(3)$, $\vartheta(4)$ etc was and if there was a clean formula for the general case. I couldn't find any general formula for $\vartheta(n)$, but WolframAlpha gave many of the desired results (such as $\vartheta(3)$ and $\vartheta(4)$). I don't know about odd $n$, but if you try out higher even $n$ you see a very apparent pattern emerging. I claim that:

$$\vartheta(2n) = (n+1)\zeta(2n+1)-\sum_{k=1}^{n-1}\zeta(2k)\zeta(2n-2k+1)$$

You can check out for yourself with e.g. $\vartheta(50)$

It is also interesting to note that:

$$\lim_{n\to\infty}\vartheta(n)=1$$

Questions

1. Is my formula of for $\vartheta(2n)$ correct?
2. How do you prove cases when $n>2$? For example $\vartheta(3) = \frac{\pi^4}{72}$?
3. Is there any other way than just numerically analyzing the answer to prove my formula?
4. Is there a formula for odd $n$ or the complete general case (for all $n$)?