If you massage your problem a little it becomes easier for Mathematica.
$f_r(x) = \sum_{n=0}^{\lfloor \frac{r-1}{2} \rfloor}c_{r,n} x^{r-2n-1} = \sum_{n=0}^{\lfloor \frac{r-1}{2} \rfloor} \frac{ (2 r - 2 n)!!}{(r- 2n-1)!} x^{r-2n-1} $.
The matrix element $X_{pq} = \int_{-1}^2 f_p(x) f_q(x) \mathrm{d} x = \sum_{n=0}^{\lfloor \frac{p-1}{2} \rfloor} \sum_{m=0}^{\lfloor \frac{q-1}{2} \rfloor} c_{p,n} c_{q,m} \int_{-1}^2 x^{p+q -2(n+m+1)} \mathrm{d} x$. After carrying out trivial integration, this becomes: $ X_{pq} = \sum_{n=0}^{\lfloor \frac{p-1}{2} \rfloor} \sum_{m=0}^{\lfloor \frac{q-1}{2} \rfloor} c_{p,n} c_{q,m} \frac{2^{p+q -2(n+m) -1} + (-1)^{p+q}}{p+q -2(n+m) -1} $
The code that would do this is as follows:
c[r_, n_] := ((2*r - 2*n)!!/(r - 2*n - 1)!); X2 = Table[ Sum[c[p, n] c[q, m] (2^(p + q - 2 (n + m) - 1) + (-1)^(p + q))/( p + q - 2 (n + m) - 1), {n, 0, Floor[(p - 1)/2]}, {m, 0, Floor[(q - 1)/2]}], {p, 5, 500}, {q, 5, 500}]
Some words of caution: The generation is slow. Generating matrix with ranges {p,5,100}, {q,5,100}
takes over 4 minutes on my machine. The matrix $X$ is manifestly symmetric, which would lead to a speed-up, reducing the run-time to about 90 seconds:
Xmat[{imin_, imax_}] := Module[{mat}, mat = PadRight[ Table[Sum[(2^(-n + p) Gamma[1 - n + p])/Gamma[-2 n + p] ( 2^(-m + q) Gamma[1 - m + q])/ Gamma[-2 m + q] ((2)^(p + q - 2 (n + m) - 1) + (-1)^(p + q))/( p + q - 2 (n + m) - 1), {n, 0, Floor[(p - 1)/2]}, {m, 0, Floor[(q - 1)/2]}], {p, imin, imax}, {q, imin, p}]]; mat + Transpose[LowerTriangularize[mat, -1]] ]