ell = 8; m = 0;
Integrate[ SphericalHarmonicY[ell, m, \[Theta], 0]* Sin[\[Theta]]^2 (1 + q Sin[\[Theta]]^2)*Sin[\[Theta]], {\[Theta], 0, \[Pi]}]

俺の環境 (ver10.4, Windows7_64bit) でこの計算やらせると落ちてしまう．
ell=6だと回る．
なんで〜〜（泣）