I was using Mathematica to get the series solution for Legendre equation. But when I get the recurrence relation and use RSolve:

RSolve[Subscript[a, n + 2] == ((n - l) (n + l + 1))/((1 + n) (2 + n)) Subscript[a, n],
Subscript[a, n], n] // FullSimplify // TraditionalForm


(where l and n are both non-negative)

Just for calculating the coefficient of the series, I only find the solution is strange:

{{Subscript[a, n]->Γ(1/2 (-l+n))/Γ(-(l/2)) (2^n (Subscript[c, 2] (-1)^n+Subscript[c, 1]) Γ(1/2 (l+n+1)))/(Γ(n+1) Γ(1/2 (l+1)))}}

I was confused that when l and n are even and n < l,the answer will be ComplexInfinity. But it is so counter-intuitive！

-
What happens after applying FullSimplify[]? – J. M. May 11 '15 at 8:00

I think you need to use this form (subscripts removed because they're messy):

FullSimplify[RSolve[a[n + 2] == ((n + k + 1) (n + k) -
l (l + 1))/((n + k + 1) (n + k + 2)) a[n], a[n], n]]

(* you should get...*)
result = (2^n (C[1]+(-1)^n C[2])Gamma[1+k]Gamma[1/2 (k-l+n)]Gamma[
1/2 (1+k+l+n)])/(Gamma[(k-l)/2] Gamma[1/2 (1+k+l)]Gamma[1+k+n])


Using this, you can now set k=0 to get the even half of the coefficients, and k=1 to get the odd half. This is the reason you were getting weird results, because you were trying to calculate the other half of coefficients from the wrong recurrence relation. NOTE: there are generated parameters C[1], C[2] because you haven't got specific initial conditions.

Now I think you can do:

evenCoeff = Table[FullSimplify[result /. {k -> 0}], {n, 0, 6, 2}]
oddCoeff = Table[FullSimplify[result /. {k -> 1}], {n, 0, 6, 2}]


Look at equations (4.13) in this text for more. I'm afraid I've never covered the Legendre equation or its polynomials so I'm only going on what I've read. The page on mathworld.wolfram.com also says this which is suggestive of you get undefined results if you try to look for the even coefficients with the odd reccurence relation and vice versa:

If l is an even integer, the series y_1(x) reduces to a polynomial of degree l with only even powers of x and the series y_2(x) diverges. If l is an odd integer, the series y_2(x) reduces to a polynomial of degree l with only odd powers of x and the series y_1(x) diverges.

-
Yes,I just think that Gamma may appear oddly when n is not necessary positive. – WateSoyan May 11 '15 at 10:18
I'm still not sure about the correctness of my answer but I think that's approximately the reason why it was giving complex infinity. Gamma returns complex infinity on all non-positive integers as well. You could try FullSimplify[... , Element[n|l,Integers]&&n>0&&l>0] see if that helps. – Histograms May 11 '15 at 10:35

Note function LegendreP gives normalized orthogonal polynomials.

Approaching your recursion by nesting can reproduce these (though not proving where coefficients vanish). So for illustrative purposes (changing l to j):

f[n_, j_] := - (j + n + 1) (j - n)/((n + 1) (n + 2));
es[a_, j_, n_] :=
First@Nest[{f[#[[2]], j] #[[1]], #[[2]] + 2} &, {a, 0}, n/2]
os[a_, j_, n_] :=
First@Nest[{f[#[[2]], j] #[[1]], #[[2]] + 2} &, {a, 1}, Floor[n/2]];
esol[x_, j_?EvenQ] := (es[1, j, #] & /@ Range[0, 2 j, 2]).Table[
x^(2 k), {k, 0, j}];
osol[x_, j_?OddQ] := (os[1, j, #] & /@ Range[1, 2 j + 1, 2]).Table[
x^(2 k + 1), {k, 0, j}];
lp[n_?OddQ, x] := With[{c = osol[1, n]}, Expand[osol[x, n]/c]]
lp[n_?EvenQ, x] := With[{c = esol[1, n]}, Expand[esol[x, n]/c]]


e.g

Grid[{Expand@LegendreP[#, x], lp[#, x]} & /@ Range[8], Frame -> All]


-