Generating function of finite sum
I would like to take the opportunity of this self answer to generalize the brilliant idea of Yurij S to calculate the generating function (g.f.) of the sum of interest. I suggest that a closed form of the g.f. of a finite (Euler) sum might be considered as a kind of substitute for the original closed form of the sum which might be lacking.
Here, as an example, we calculate directly the generating function of the original sum
$$p_2(n)=\sum_{k=1}^n \frac{H_k}{2k+1}$$
hoping to get a closed expression for the g.f. (the hope turns out to be justified).
The g.f. of our sum of interest is defined as
$$g_2(z) = \sum_{k=1}^\infty p_2(n) z^n$$
We proceed in four steps.
Step 1: write sum as double integral
Using the relations $H_k = \int_0^1 \frac{1-x^k}{1-x}\,dx$ and $\frac{1}{2k+1}=\int_0^1 y^{2k}\,dy$ the sum $p_2(n)$ transforms itself naturally into a double integral:
$$p_2(n) = \int_0^1 \,dx \int_0^1 \,dy \sum_{k=1}^n y^{2k} \frac{1-x^k}{1-x}$$
Step 2: perform the finite sum
Doing the sum under the integral gives the integrand
$$i(n) = \frac{y^2 \left(\left(y^2-1\right) x^{n+1} y^{2 n}-x y^{2 n+2}+y^{2 n}+x-1\right)}{(x-1) \left(y^2-1\right) \left(x y^2-1\right)}$$
Step 3: form the g.f. with the Integrand
Before integrating we first form the g.f. under the integral
$$i_g(z) = \sum_{n=1}^\infty i(n) z^n = -\frac{y^2 z}{(z-1) \left(y^2 z-1\right) \left(x y^2 z-1\right)}$$
Step 4: do the integration (in an appropriate order)
Now we do the $x$-integral
$$i_{g,x}(z)= \int_0^1 i_g(z) \,dx = -\frac{\log \left(1-y^2 z\right)}{(1-z) \left(1-y^2 z\right)}$$
and finally the $y$-integral gives the g.f.
$$g(z)= \int_0^1 i_{g,x}(z) \,dy = \frac{1}{12 (z-1) \sqrt{z}} \left(-12 \text{Li}_2\left(\frac{\sqrt{z}-1}{\sqrt{z}+1}\right)+12 \tanh ^{-1}\left(\sqrt{z}\right)^2+12 \left(\log (4-4 z)-2 \log \left(\sqrt{z}+1\right)\right) \tanh ^{-1}\left(\sqrt{z}\right)-\pi ^2\right)$$
This is the desired closed form of the g.f.
The series expansion of $g(z)$ about $z=0$ starts like this
$$g(z) = \frac{z}{3} + \frac{19 z^2}{30}+ \frac{94 z^3}{105}+\frac{4259 z^4}{3780}+\frac{2774 z^5}{2079}+...$$
It has only positive integer powers of $z$ as it should, and it is easy to identify the coefficient of $z^n$ as the value of $p_2(n)$ by direct comparison (left as an easy exercise to the reader).
Lemma: relation between generating functions
There is an interesting and useful relation between the g.f. of a series $a(n)$ and the g.f. of a finite sum of the series.
Let
$$g_a(z) = \sum_{n=1}^\infty a(n) z^n$$
$$s(n) = \sum_{k=1}^n a(k)$$
$$g_s(z) = \sum_{n=1}^\infty s(n) z^n$$
then
$$g_s(z) = \frac{1}{1-z} g_a(z)$$
The proof is easily found by interchanging the order of summation in $g_s(z)$
Here are some simple applications of the lemma.
Harmonic number
Taking $a(k)=1/k$ gives the generating function $g_a(z) = \sum_{k=1}^\infty \frac{z^k}{k} = - \log(1-z)$. Hence for the sum we get the g.f. $g_s(z) = \frac{g_a(z)}{1-z} = \frac{ - \log(1-z)}{1-z}$. Noticing thet the sum is just the harmonic number we have recovered the g.f. of the latter.
Iterated sums
Defining the $q$-fold iterated sum recursively by
$s(0,n) = a(n)$, $s(q+1,n) = \sum_{k=1}^n s(q,k)$
the g.f. of $s(q,n)$ is given by the simple formula
$$g(q,z) = \sum_{k=1}^\infty z^n s(q,n) = \frac{g(0,z)}{(1-z)^q}$$