where $H_k$ is the $k$-th harmonic number, $H_k=\sum_{j=1}^k H_k=\sum\limits_{j=1}^k \frac1{j}$.
Might there be a (hopefully simple) $n$-term recurrence that generates these polynomials?
Addendum:
After staring long and hard at Pietro's answer, I feel now that my second question was sorta kinda dumb; I already knew that both Legendre functions satisfied the same difference equation, so it stands to reason that a linear combination of them should also be a solution to that recurrence.
I now would like to expand my first question a bit: is it possible to express $W_n(z)$ as a single hypergeometric function (e.g. ${}_p F_q$ or some of the fancy multivariate ones), perhaps with one of the parameters being a negative integer? For instance, $P_n(z)$ is expressible as a Gaussian hypergeometric function ${}_2 F_1$ with one of the numerator parameters being a negative integer. Might there be something similar for the $W_n$?
I would also like to consider an additional question: are the $W_n$ orthogonal polynomials with respect to some weight function $\omega(x)$ and an associated support interval $(a,b)$? That is, if
$$\int_a^b\omega(t)W_j(t)W_k(t)\mathrm dt=0,\qquad j\neq k$$
for some $\omega(x)$ and some interval $(a,b)$, what is this weight function and its support interval?

