I encountered the following integral in my (physics) research, and I've yet to find an analytic solution: $I(n_1,n_2,n_3) = \int_{-1}^{1} d(\cos\theta_1) \int_{-1}^{1} d(\cos\theta_2) P_{n_1}(\cos\theta_1) P_{n_2}[\cos(\theta_1-\theta_2)] P_{n_3}(\cos\theta_2)$ where $P_n(x)$ is the nth Legendre polynomial. For reference, the orthogonal polynomial relation for Legendre polynomials is: $\int_{-1}^{1} d(\cos\theta) P_{n_1}(\cos\theta) P_{n_2}(\cos\theta) = \int_{0}^{\pi} d\theta (\sin\theta) P_{n_1}(\cos\theta) P_{n_2}(\cos\theta) = \frac{2}{2n_1+1}\delta_{n_1,n_2}$
Mathematica does not have much trouble in evaluating the integral exactly (Integrate, not NIntegrate) for a given $\{n_1, n_2, n_3\}$.
k[n1_,n2_,n3_]:=Integrate[Sin[x]Sin[y]LegendreP[n1,Cos[x]]LegendreP[n2,Cos[x-y]]LegendreP[n3,Cos[y]],{x,0,Pi},{y,0,Pi}]
I can see already that $I(n_1, n_2, n_3)$ is not diagonal. Even the diagonal terms contain factors which I've yet to figure out, here's what I see so far: $I(n,n,n) = 2^n \left[\frac{2}{2n+1}\right]^2 \times \left\{1,\frac{1}{2},\frac{1}{3},\frac{1}{5},\frac {4}{35},\frac{4}{63},\frac{8}{231},\frac{8}{429}, \frac{64}{6435}, ... \right\}$ where the bracketed terms correspond to $n=0,1,2,...$
Many of the off-diagonal terms appear to contain a $\pi^2$ term and some power of 2.
Do you have any suggestions? Useful variable transformations, pattern recognition? Thanks in advance.