A recurrence relation for the Harmonic numbers of the form $H_n = \sum\limits_{k=1}^{n-1}f(k,n)H_k$

It helps to visualize the terms in a square of products $1/(ij)$ with $i$ and $j$ running for $1$ to $n$. The sum over the first term contains all products with $i\ne j$ exactly once, whereas the sum over the second term roughly corresponds to the upper left half of the square, but with the left-most column, which adds up to $H_n$, excluded. Thus we have

$$ \def\sub#1{{\scriptstyle{i\ne j}\atop{\scriptstyle i,j\le #1}}} \sum_{k=1}^{n-1}\frac{2}{k+1}H_k=\sum_{\sub n}\frac1{ij}$$

and

$$-\sum_{k=1}^{n-1}\frac{1}{1+n-k}=H_n-\sum_{i+j\le n+1}\frac1{ij}\;.$$

Substituting this into your equation, multiplying through by $n-1$ and simplifying leads to

$$\sum_{i+j\le n+1}\frac1{ij}-\sum_{\sub n}\frac1{ij}=\frac2{n+1}H_n\;,$$

$$\sum_{i+j\le n+1}\frac1{ij}-\sum_{\sub n}\frac1{ij}=2\sum_i\frac1{n+1}\frac1i\;,$$

$$\sum_{i+j\le n+1}\frac1{ij}=\sum_{\sub{n+1}}\frac1{ij}\;.$$

This we can prove by induction: The equation is satisfied for $n=0$, and going from $n$ to $n+1$ adds

$$\sum_{i+j=n+1}\frac1{ij}=\sum_{i=1}^n\frac1i\frac1{n+1-i}$$

to the left-hand side and also

$$2\frac1{n+1}\sum_{i=1}^n\frac1i=\frac1{n+1}\sum_{i=1}^n\left(\frac1i+\frac1{n+1-i}\right)=\sum_{i=1}^n\frac1i\frac1{n+1-i}$$

to the right-hand side.

Here's a generating function route: as I already mentioned in the comments,

$$\begin{align*} H_n&=\frac{n+1}{n-1}\sum_{k=1}^{n-1}\left(\frac{2}{k+1}-\frac{1}{n-k+1}\right)H_k\\ H_n&=\frac{n+1}{n-1}\left(2\sum_{k=1}^{n-1}\frac{H_k}{k+1}-\sum_{k=1}^{n-1}\frac{H_k}{n-k+1}\right)\\ H_n&=\frac{n+1}{n-1}\left(2\sum_{k=1}^{n}\frac{H_k}{k+1}-\sum_{k=1}^{n}\frac{H_k}{n-k+1}+H_n-\frac{2H_n}{n+1}\right)\\ H_n&=\frac{n+1}{n-1}\left(2\sum_{k=1}^{n}\frac{H_k}{k+1}-\sum_{k=1}^{n}\frac{H_k}{n-k+1}\right)+\frac{n+1}{n-1}H_n-\frac2{n-1}H_n\\ H_n&=\frac{n+1}{n-1}\left(2\sum_{k=1}^{n}\frac{H_k}{k+1}-\sum_{k=1}^{n}\frac{H_k}{n-k+1}\right)+H_n\\ \sum_{k=1}^{n}\frac{H_k}{n-k+1}&=2\sum_{k=1}^{n}\frac{H_k}{k+1} \end{align*}$$

We note that the sum $\sum\limits_{k=1}^{n}\frac{H_k}{n-k+1}$ is in the form of a convolution; thus, its generating function is

$$\left(\sum_{j=1}^\infty H_j x^j\right)\left(\sum_{j=1}^\infty \frac{x^j}{j}\right)=\frac{(\log(1-x))^2}{1-x}$$

The remaining task is to prove that the generating function given above is also the generating function of $2\sum\limits_{k=1}^{n}\frac{H_k}{k+1}$; to that effect, there is the identity

$$2\sum_{k=1}^{n}\frac{H_k}{k+1}=(H_{n+1})^2-H_{n+1}^{(2)}$$

where $H_n^{(k)}=\sum\limits_{j=1}^n \frac1{j^k}$ is a generalized harmonic number. From here and here (see formula 36), we have the generating functions

$$\begin{align*} \sum_{j=1}^\infty H_{j+1}^{(2)}x^{j+1}&=\frac{\mathrm{Li}_2(x)}{1-x}-x\\ \sum_{j=1}^\infty (H_{j+1})^2 x^{j+1}&=\frac{(\log(1-x))^2+\mathrm{Li}_2(x)}{1-x}-x \end{align*}$$

where $\mathrm{Li}_2(x)=-\int_0^x \frac{\log(1-u)}{u}\mathrm du$ is a dilogarithm.

Thus,

$$\sum_{j=1}^\infty (H_{j+1})^2 x^{j+1}-\sum_{j=1}^\infty H_{j+1}^{(2)}x^{j+1}=\frac{(\log(1-x))^2}{1-x}$$

and we're golden.

Thanks for your help. Here is another proof that I found today :

From the well known $\zeta(3)=\frac{1}{2}\sum_{k=1}^{\infty}\frac{H_k}{k^2}$, we find $$ \zeta(3)=\sum_{k=1}^{\infty}\frac{H_k}{(k+1)^2} \quad\quad(1) $$ Were $\zeta(z)$ is the Riemann zeta function. But also, $$ \zeta(3)=\frac{1}{2}\int_{0}^{\infty}\frac{t^2}{e^t-1}dt=4\int_{0}^{\pi/2}\tan{x}(\ln{\sin{x}})^2dx=4\int_{0}^{\pi/2}\tan{x}\left(\sum_{k=1}^{\infty}\frac{\cos^{2k}{x}}{2k}\right)^2dx $$ Where I used the substitution $e^{-t}=\sin^2{x}$. By rearranging and using the formula for raising power series to powers (e.g. 0.314 p.17 in Gradshteyn and Ryzhik's Table of Integrals, Series, and Products), $$ \zeta(3)=\sum_{k=0}^{\infty}c_k \int_{0}^{\pi/2}\sin{x}\ \cos^{2k+3}{x}dx=\sum_{k=0}^{\infty}c_k \frac{1}{2}B(1,k+2)=\sum_{k=1}^{\infty}\frac{c_{k-1}}{2(k+1)} $$ $$ \text{where,}\quad c_0=1,\quad c_k=\frac{1}{k}\sum_{i=1}^{k}\frac{3i-k}{i+1}c_{k-i} $$ Equating with (1) and rearranging, we get desired relation.