A very long comment, definitely not a "simple real method" that OP is looking for.

The structure of such sums in terms of multiple zeta values (MZVs) are considered well-understood. Several years ago I implemented an algorithm to calculate logarithmic integrals via MZVs, the algorithm for Euler sums is almost the same, but I procrastinated on implementing them until recently. Now you can compute them using the same Mathematica package.

An overview of the theory and algorithm behind can be learned here.

For OP's sum, one simply input the program only accepts sum from $0$ to $\infty$, so we need to do a shift by $1$. This corroborates OP's result.

What OP called "skew-harmonic number" $\overline{H}_n = \sum_{i=1}^n \frac{(-1)^{i-1}}{i}$ is represented in the program as MultiHarmonic[{1},{1},{-1},n], the two sums OP mentioned can also be computed using same algorithm:

OP asks about the sum $\sum_{n=1}^{\infty}\frac{H_n H_{4n}}{(2n-1)^2}$ at the end, it equals note that the weight is now inhomogeneous, in contrast to all previous examples.

All variations of Euler sums in the literature can be routinely converted to MZVs. More examples that OP might enjoy solving using "simple real method" are: $$\sum_{n=1}^\infty \frac{(H_n)^7}{n^2} \qquad \sum_{n=1}^\infty \frac{H_n H_n^{(2)} H_n^{(3)}}{n^3}\qquad \sum_{n=1}^\infty \frac{(-1)^n \overline{H}_n \overline{H}_n^{(2)} H_n^{(2)}}{n}$$

$$\sum_{n=0}^\infty \frac{H_n H_{2n}}{(4n+1)^2} \qquad \sum_{n=0}^\infty \frac{H_n H_{4n}}{(4n+1)^2}\qquad \sum_{n=0}^\infty \frac{H_{2n} H_{4n}}{(4n+1)^2}$$ the middle sum, for example, is $$-\frac{\pi ^2 G}{32}+\frac{7}{8} G \log ^2(2)+\frac{\pi}{8} G \log (2)-\beta(4)+\frac{\pi}{4} \Im \text{Li}_3(\frac{1+i}{2})+\Im\text{Li}_4(\frac{1+i}{2})+2 \log (2) \Im \text{Li}_3(\frac{1+i}{2})+\frac{3 \text{Li}_4\left(\frac{1}{2}\right)}{2}+\frac{7 \pi \zeta (3)}{64}+\frac{21}{32} \zeta (3) \log (2)-\frac{53 \pi ^4}{3840}+\frac{\log ^4(2)}{16}-\frac{11}{192} \pi \log ^3(2)-\frac{9}{128} \pi ^2 \log ^2(2)-\frac{9}{256} \pi ^3 \log (2)$$

as well as $$\sum_{n\geq 0} \frac{(-1)^n A_n H_n^{(2)} H_n}{2n+1} \qquad \sum_{n\geq 0} \frac{(-1)^n A_n H_n^{(2)}}{n+1} \qquad A_n = \sum_{i=0}^{n-1} \frac{(-1)^{i-1}}{2i-1}$$

Similar to logarithmic integrals, one can literally continue forever: there are millions if not billions of variations, they can however all be tackled using the same method.

Nonetheless, there are still much room of improvement for the program: currently it can't

1. reduce sums such as $\sum_{n\geq 1} \frac{H_n^{(p)}}{n^q}$ with $p+q$ odd using $\zeta(s)$, this is well-known to be possible (parity-depth reduction of MZV)
2. handle generalizations of Euler sums that contain binomial coefficient such as this one.

Algorithms for both points are only recently understood, and will be implemented in future.

A very long comment, definitely not a "simple real method" that OP is looking for.

The structure of such sums in terms of multiple zeta values (MZVs) are considered well-understood. Several years ago I implemented an algorithm to calculate logarithmic integrals via MZVs, the algorithm for Euler sums is almost the same. Now you can compute them using the same Mathematica package.

An overview of the theory and algorithm behind can be learned here.

For OP's sum, one simply input the program only accepts sum from $0$ to $\infty$, so we need to do a shift by $1$. This corroborates OP's result.

What OP called "skew-harmonic number" $\overline{H}_n = \sum_{i=1}^n \frac{(-1)^{i-1}}{i}$ is represented in the program as MultiHarmonic[{1},{1},{-1},n], the two sums OP mentioned can also be computed using same algorithm:

OP asks about the sum $\sum_{n=1}^{\infty}\frac{H_n H_{4n}}{(2n-1)^2}$ at the end, it equals note that the weight is now inhomogeneous.

All variations of Euler sums in the literature can be routinely converted to MZVs. More examples that OP might enjoy solving using "simple real method" are: $$\sum_{n=1}^\infty \frac{(H_n)^7}{n^2} \qquad \sum_{n=1}^\infty \frac{H_n H_n^{(2)} H_n^{(3)}}{n^3}\qquad \sum_{n=1}^\infty \frac{(-1)^n \overline{H}_n \overline{H}_n^{(2)} H_n^{(2)}}{n}$$

$$\sum_{n=0}^\infty \frac{H_n H_{2n}}{(4n+1)^2} \qquad \sum_{n=0}^\infty \frac{H_n H_{4n}}{(4n+1)^2}\qquad \sum_{n=0}^\infty \frac{H_{2n} H_{4n}}{(4n+1)^2}$$ the middle sum, for example, is $$-\frac{\pi ^2 G}{32}+\frac{7}{8} G \log ^2(2)+\frac{\pi}{8} G \log (2)-\beta(4)+\frac{\pi}{4} \Im \text{Li}_3(\frac{1+i}{2})+\Im\text{Li}_4(\frac{1+i}{2})+2 \log (2) \Im \text{Li}_3(\frac{1+i}{2})+\frac{3 \text{Li}_4\left(\frac{1}{2}\right)}{2}+\frac{7 \pi \zeta (3)}{64}+\frac{21}{32} \zeta (3) \log (2)-\frac{53 \pi ^4}{3840}+\frac{\log ^4(2)}{16}-\frac{11}{192} \pi \log ^3(2)-\frac{9}{128} \pi ^2 \log ^2(2)-\frac{9}{256} \pi ^3 \log (2)$$

as well as $$\sum_{n\geq 0} \frac{(-1)^n A_n H_n^{(2)} H_n}{2n+1} \qquad \sum_{n\geq 0} \frac{(-1)^n A_n H_n^{(2)}}{n+1} \qquad A_n = \sum_{i=0}^{n-1} \frac{(-1)^{i-1}}{2i-1}$$

Similar to logarithmic integrals, one can literally continue forever: there are millions if not billions of variations, they can however all be tackled using the same method.

Nonetheless, there are still much room of improvement for the program: currently it can't

1. reduce sums such as $\sum_{n\geq 1} \frac{H_n^{(p)}}{n^q}$ with $p+q$ odd using $\zeta(s)$, this is well-known to be possible (parity-depth reduction of MZV)
2. handle generalizations of Euler sums that contain binomial coefficient such as this one.

Algorithms for both points are understood, and will be implemented in future.