Please don't be frightened by the length of this, I just wanted to provide ample detail. If you want, you can skip the derivation and go straight to the result at the bottom.
I have
$$ \begin{multline}\int_{-1}^{1}\left(\sum_{n_1-l_1=0}^{\infty}C_{n_1-l_1}^{l_1+1}(x)r^{n_1-l_1}\right)\left(\sum_{n_2-l_2=0}^{\infty}C_{n_2-l_2}^{l_2+1}(x)r^{n_1-l_1}\right)\times\\ \left(\sum_{n_3-l_3=0}^{\infty}C_{n_3-l_3}^{l_3+1}(x)r^{n_3-l_3}\right)(1-x^2)^{(l_1+l_2+l_3+1)/2}dx \end{multline} $$ Applying the identity $\left(\sum_{n=0}^\infty a_nr^n\right)\left(\sum_{n=0}^\infty b_nr^n\right)=\sum_{n=0}^\infty\sum_{k=0}^n a_kb_{n-k}r^n$ twice to the above gives: $$ =\sum_{i=0}^{\infty}\sum_{j=0}^{i}\sum_{k=0}^{i-j}r^{i}\int_{-1}^{1}C_{j}^{l_1+1}(x)C_{k}^{l_2+1}(x)C_{i-j-k}^{l_3+1}(x)(1-x^2)^{(l_1+l_2+l_3+1)/2}dx $$
Using the generating function for the Gegenbauer Polynomials: $$ (1-2xr+r^2)^{-\lambda}=\sum_{n=0}^{\infty}C_{n}^{\lambda}(x)r^n\ \ \ ,\left|r\right|<1 $$ And writing $L=l_1+l_2+l_3$, we get: $$ =\int_{-1}^{1}(1-2xr+r^2)^{-(L+3)}(1-x^2)^{(L+3)/2} $$ $$ =\sqrt{\pi}\frac{\Gamma(\frac{L+3}{2})}{\Gamma(\frac{L+4}{2})}(1-r^2)^{-(L+3)} $$ Rewriting $(1-r^2)^{-(L+3)}=(1-r^2)^{-(l_1+1)}(1-r^2)^{-(l_2+1)}(1-r^2)^{-(l_3+1)}$ and applying the binomial theorem gives: $$ \begin{multline} \sqrt{\pi}\frac{\Gamma\left(\frac{L+3}{2}\right)}{\Gamma\left(\frac{L+4}{2}\right)}\left(\sum_{n=0}^{\infty}\left(-1\right)^n\left(\begin{array}{c}-l_1-1\\n\end{array}\right)r^{2n}\right)\times\\ \left(\sum_{n=0}^{\infty}\left(-1\right)^n\left(\begin{array}{c}-l_2-1\\n\end{array}\right)r^{2n}\right)\left(\sum_{n=0}^{\infty}\left(-1\right)^n\left(\begin{array}{c}-l_3-1\\n\end{array}\right)r^{2n}\right) \end{multline} $$ Applying the product of series identity twice again gives: $$ =\sum_{n=0}^{\infty}\sum_{l=0}^{n}\sum_{m=0}^{n-l}r^{2n}(-1)^n\left(\begin{array}{c}-l_1-1\\l\end{array}\right)\left(\begin{array}{c}-l_2-1\\m\end{array}\right)\left(\begin{array}{c}-l_3-1\\n-l-m\end{array}\right)\sqrt{\pi}\frac{\Gamma\left(\frac{L+3}{2}\right)}{\Gamma\left(\frac{L+4}{2}\right)} $$
Comparing the series in $i,j,k$ to this one in $n,l,m$ I find that $i=2n$ (so that the original series is zero unless i is an even integer). Then comparing the coefficients of $r$ I get: $$ \begin{multline} \sum_{j=0}^{i}\sum_{k=0}^{i-j}\int_{-1}^{1}C_{j}^{l_1+1}(x)C_{k}^{l_2+1}(x)C_{i-j-k}^{l_3+1}(x)(1-x^2)^{(l_1+l_2+l_3+1)/2}dx\\=\sum_{l=0}^{n}\sum_{m=0}^{n-l}(-1)^n\left(\begin{array}{c}-l_1-1\\l\end{array}\right)\left(\begin{array}{c}-l_2-1\\m\end{array}\right)\left(\begin{array}{c}-l_3-1\\n-l-m\end{array}\right)\sqrt{\pi}\frac{\Gamma\left(\frac{L+3}{2}\right)}{\Gamma\left(\frac{L+4}{2}\right)} \end{multline} $$ I then set $$ \begin{align} i=2n\\j=l=n_1-l_1\\k=m=n_2-l_2\\i-j-k=n_3-l_3 \end{align} $$ Whence I get the $\bf{FINAL\ RESULT}$: $$ \int_{-1}^{1}C_{n_1-l_1}^{l_1+1}(x)C_{n_2-l_2}^{l_2+1}(x)C_{n_3-l_3}^{l_3+1}(x)(1-x^2)^{(l_1+l_2+l_3+1)/2}dx= $$ $$ (-1)^{(N-L)/2}\left(\begin{array}{c}-l_1-1\\ n_1-l_1\end{array}\right)\left(\begin{array}{c}-l_2-1\\n_2-l_2\end{array}\right)\left(\begin{array}{c}-l_3-1\\-\frac{1}{2}((n_1-l_1)+(n_2-l_2)-(n_3-l_3))\end{array}\right)\sqrt{\pi}\frac{\Gamma\left(\frac{L+3}{2}\right)}{\Gamma\left(\frac{L+4}{2}\right)} $$
So, $\bf{HERE's}$ my question:
Why is this incorrect? I've tested it and it is not true for all values of $n_1,l_1,n_2,l_2,n_3,l_3$.
I've now spent the better part of a week trying to figure out why this doesn't work, but I can't figure it out. It seems like maybe it has something to do with the way I equate the indices, or perhaps a theorem about the convergence of the product of power series that would make something I did illegal, but I can't find anything so far that would tell me that there's anything wrong with it. In case it's helpful the answer seems to be off by a multiplicative factor that is a function of the indices (a gamma function maybe?).