I might as well: there is an algorithm, due to Salzer, for converting to and from different orthogonal polynomial expansions. (I previously talked about this algorithm here.) If you intend to do this entire business of conversions, I believe an algorithm might be more useful to you than an integral expression.
What follows are a few relevant Mathematica routines. Even though I used Mathematica, I hope that the algorithm is still transparent and easily translatable (the indexing starts from 1; some adjustment is necessary if the indexing in your language starts at 0):
monomialToLegendre[pc_?VectorQ] :=
Module[{n = Length[pc] - 1, lc, q},
lc[1] = pc[[n]]; lc[2] = Last[pc];
Do[
q[1] = lc[1];
lc[1] = pc[[n - k + 1]] + lc[2]/3;
Do[
q[j] = lc[j];
lc[j] = j lc[j + 1]/(2 j + 1) + (j - 1) q[j - 1]/(2 j - 3);
, {j, 2, k - 1}];
q[k] = lc[k];
lc[k] = (k - 1) q[k - 1]/(2 k - 3);
lc[k + 1] = k q[k]/(2 k - 1);
, {k, 2, n}];
Table[lc[k], {k, n + 1}]]
One could use Salzer again for converting from the Legendre to the monomial basis. However, it is a bit clearer (and gives rise to essentially the same method) to treat the problem as one of Taylor/Maclaurin expansion of the polynomial expressed in terms of Legendre polynomials. One can then use Clenshaw's algorithm for the purpose:
legendreToMonomial[lc_?VectorQ] :=
Module[{n = Length[lc] - 1, pc, v = 0, w, z},
pc[1] = Last[lc]; pc[_] = z[_] = 0;
Do[
w = pc[1];
pc[1] = lc[[j]] - v z[1];
z[1] = w;
Do[
w = pc[i];
pc[i] = (2 j - 1) z[i - 1]/j - v z[i];
z[i] = w;
, {i, 2, n - j + 2}];
v = (j - 1)/j;
, {j, n, 1, -1}];
Table[pc[j], {j, n + 1}]]
As an example, to demonstrate the identity
$$5P_0(x)-4P_1(x)-2P_2(x)+3P_3(x)+P_4(x)=\frac{35}{8}x^4+\frac{15}{2}x^3-\frac{27}{4}x^2-\frac{17}{2}x+\frac{51}{8}$$
monomialToLegendre[legendreToMonomial[{5, -4, -2, 3, 1}]] ==
{5, -4, -2, 3, 1}
True
legendreToMonomial[monomialToLegendre[{51/8, -17/2, -27/4, 15/2, 35/8}]]
== {51/8, -17/2, -27/4, 15/2, 35/8}
True
Best Answer
The most general solution of Legendre equation is $$y = A{P_n} + B{Q_n}.$$ Let $y(x) = A(x){P_n}(x)$. Then $y' = AP' + A'P$ and $y'' = AP'' + 2A'P' + A''P$. So $$(1 - {x^2})(AP'' + 2A'P' + A''P) - 2x(AP' + A'P) + n(n + 1)AP = 0.$$ Note that $$(1 - {x^2})(AP'') - 2x(AP') + n(n + 1)AP = 0$$ which means some terms in the above equation vanish. Now let $A' = u$ and reduce the order so that $$2\frac{{dP}}{P} + \frac{{du}}{u} - \frac{{2xdx}}{{1 - {x^2}}} = 0$$ and $$u = \frac{{{\text{const}}}}{{(1 - {x^2}){P^2}}}$$ so $$A = {C_n}\int {\frac{1}{{(1 - {x^2}){P^2}}}dx}.$$ See if you can do the rest.