Legendre Polynomials Triple Product I have to solve the following integral:
\begin{align}
\int_{-1}^{1} \left(x^2 -1\right)^3 P_k(x)\,P_l(x)\, P_m(x) \;dx
\end{align}
where $P_{k,l,m}$ are Legendre Polynomials
The triple product
\begin{align}
\int_{-1}^{1} P_k(x)\,P_l(x)\, P_m(x) \;dx = 2 \begin{pmatrix} k & l & m \\ 0 & 0 & 0 \end{pmatrix}^2
\end{align}
using the special case of $3j$ symbol form
\begin{align}
\begin{pmatrix} k & l & m \\ 0 & 0 & 0 \end{pmatrix}
    &= (-1)^s \sqrt{(2s-2k)! (2s-2l)! (2s-2m)! \over (2s+1)!}
        {s! \over (s-k)! (s-l)! (s-m)!}
        \\ & \mbox{for $2s=k+l+m$ even} \\[3pt]
\begin{pmatrix} k & l & m \\ 0 & 0 & 0 \end{pmatrix}
    &= 0
        \quad\mbox{for $2s=k+l+m$ odd} \\
\end{align} 
I'm sure you should be able to solve this by doing integration by parts but can't seem to get it to work. Any tips?
So using the answer below I think you get the following for step 1 of 3
\begin{multline}
 \int_{-1}^{1}(x^2-1)^3 P_k P_l P_m 
 = \overbrace{(x^2-1)^3\frac{(P_{k+1} - P_{k-1})}{2k+1} P_l P_m \Big]_{-1}^1}^\text{ = 0}\\
-\int_{-1}^{1} \frac{(P_{k+1} - P_{k-1})}{2k+1}(x^2-1)^2\Big( 
6xP_l P_m \\
+  (1+l) P_m(P_{l+1} - P_{l-1})  + 
(1+m) P_l(P_{m+1} - P_{m-1}) \Big) \; dx\\ 
\end{multline}
Not sure if the formula for integration works as I think the $6xP_lP_m$ term might cause problems?

I think the best way to approach this is as follows,
note that 
\begin{align}
(x^2 -1 ) = \frac{P_2 - 2}{3} 
\end{align}
You can then use the following definition
\begin{align}
P_kP_l = \sum_{m=|k-l|}^{k+l} 
\begin{pmatrix}
k & l & m \\
0 & 0 & 0
\end{pmatrix}^2
(2m+1)P_m
\end{align}
This allows the integral to be written as follows
\begin{align}
\int_{-1}^{1} (x^2-1)^3P_iP_jP_k \; dx &=  \int_{-1}^{1} \frac{1}{9}\left(P_2^3 + . . .-8 \right) P_i P_j P_k \; dx
\end{align}
The most difficult term to deal with is the $ P_2^3 P_i P_j P_k$ 
\begin{align}
P_2^3 P_i P_j P_k &= \sum_{m=0}^{4}
\begin{pmatrix}
2 & 2 & m \\
0 & 0 & 0
\end{pmatrix}^2 (2m+1)P_m P_2 P_i P_j P_k \\
&= \sum_{m=0}^{4}
\begin{pmatrix}
2 & 2 & m \\
0 & 0 & 0
\end{pmatrix}^2 (2m+1)
\sum_{n=|m-2|}^{m+2}
\begin{pmatrix}
2 & m & n \\
0 & 0 & 0
\end{pmatrix}^2 (2n+1)P_n P_i P_j P_k \\
&= 
\sum_{m=0}^{4}
\begin{pmatrix}
2 & 2 & m \\
0 & 0 & 0
\end{pmatrix}^2 (2m+1)
\sum_{n=|m-2|}^{m+2}
\begin{pmatrix}
2 & m & n \\
0 & 0 & 0
\end{pmatrix}^2 (2n+1)
\sum_{l=|n-i|}^{n+i}
\begin{pmatrix}
n & i & l \\
0 & 0 & 0
\end{pmatrix}^2 (2l+1)
P_l P_j P_k
\end{align}
Which can then make use of the usual triple integral.
All other terms can be solved for in a similar manner.
