What is the relationship between hypergeometric function and Legendre polynomials? I have an equation 
$$
-\frac{1}{2}y''(x)-b~\text{sech}^2(ax)~y=-\frac{k^2}{2}y
$$
and I know that it has solution in terms of Legendre polynomials:
$$
y_1(x)=P_s^\epsilon(\zeta)
$$
$$
y_2(x)=Q_s^\epsilon(\zeta)
$$
where 
$$
s=\frac{1}{2} \left(\sqrt{\frac{8 b}{a^2}+1}-1\right)
$$
$$
\epsilon =\frac{k}{a}
$$
$$
\zeta =\tanh (a x)
$$
I also know that the solution can be written in the form of Hypergeometric function:
$$
y_3=\left(1-\zeta ^2\right)^{\epsilon /2} \, _2F_1\left(\epsilon -s,s+\epsilon +1;\epsilon +1;\frac{1-\zeta }{2}\right)
$$
Since there are are only two linear independent solutions, $y_3$ must can be expressed in terms of $y_1$ and $y_2$, so what is the relationship between them? And what is the fourth solution $y_4$ which is linear independent with respect to $y_3$? 

Here are some mathematica code to verify the solution:
eqs = -(1/2) y''[x] - b (Sech[a x])^2 y[x] == -(k^2/2) y[x];
ζ = Tanh[a x]; s = 1/2 (Sqrt[1 + (8 b)/a^2] - 1); ϵ = k/a;
y1 = LegendreP[s , ϵ, ζ];
y2 = LegendreQ[s , ϵ, ζ];
y3 = (1 - ζ^2)^(ϵ/2)Hypergeometric2F1[ϵ - s, ϵ + s + 1, ϵ + 1, 1/2 (1 - ζ)];

eqs /. y -> Function[{x}, Evaluate[y1]] // FullSimplify
eqs /. y -> Function[{x}, Evaluate[y2]] // FullSimplify
eqs /. y -> Function[{x}, Evaluate[y3]] // FullSimplify

(*True*)
(*True*)
(*True*)


Indeed you can express associated Legendre functions in terms of hypergeometric functions. See, for example, Abramowitz / Stegun, section 8.1, or the Digital Library of Mathematical Functions, section 14.3:
http://people.math.sfu.ca/~cbm/aands/page_332.htm
http://people.math.sfu.ca/~cbm/aands/page_333.htm
http://dlmf.nist.gov/14.3
In your particular case, equations 14.3.15 and 14.3.19 of the latter should give you the answer:
http://dlmf.nist.gov/14.3#E15
http://dlmf.nist.gov/14.3#E19
