Efficient evaluation of a hyper-spherical harmonic expansion.

I need to integrate a function that is expressed as a hyper-spherical harmonic expansion. This function has many terms and the basis functions themselves involve many operations, as such evaluation of the function is prohibitively slow for integration. As an example here is one of the (smaller) basis functions (of which there are many, and yes, I've tried to simplify the basis functions as much as possible with a computer algebra package):

(0.44946657.*(154.^(-1/2).*pi.^(-1).*(1+3.*cos(2.*t)).*(50+90.*cos(w)+63.*cos(2.*w)+28.*cos(3.*w)).*sin((1/2).*w).^2)+0.04407382.*((3/8).*1001.^(-1/2).*pi.^(-1).*(9+20.*cos(2.*t)+35.*cos(4.*t)).*(45+70.*cos(w)+28.*cos(2.*w)).*sin((1/2).*w).^4)-0.05214885.*(3.*(5/143).^(1/2).*pi.^(-1).*cos(4.*p).*(45+70.*cos(w)+28.*cos(2.*w)).*sin(t).^4.*sin((1/2).*w).^4)+0.55470020.*(6.*(70/13).^(1/2).*pi.^(-1).*cos(t).*cos((1/2).*w).*(5+8.*cos(w)).*sin(4.*p).*sin(t).^4.*sin((1/2).*w).^5)+0.26590801.*((1/16).*110.^(-1/2).*pi.^(-1).*(50+105.*cos(2.*t)+126.*cos(4.*t)+231.*cos(6.*t)).*(7+8.*cos(w)).*sin((1/2).*w).^6)+0.10050378.*(3.*(7/110).^(1/2).*pi.^(-1).*cos(4.*p).*(9+11.*cos(2.*t)).*(7+8.*cos(w)).*sin(t).^4.*sin((1/2).*w).^6)-0.43852901.*(12.*(7/13).^(1/2).*pi.^(-1).*cos(t).*(7+13.*cos(2.*t)).*cos((1/2).*w).*sin(4.*p).*sin(t).^4.*sin((1/2).*w).^7)+0.29985590.*(715.^(-1/2).*pi.^(-1).*(35+(-1260).*cos(t).^2+6930.*cos(t).^4+(-12012).*cos(t).^6+6435.*cos(t).^8).*sin((1/2).*w).^8)-0.03268602.*((3/4).*(7/65).^(1/2).*pi.^(-1).*cos(4.*p).*(99+156.*cos(2.*t)+65.*cos(4.*t)).*sin(t).^4.*sin((1/2).*w).^8)-0.34860834.*(3.*pi.^(-1).*cos(8.*p).*sin(t).^8.*sin((1/2).*w).^8))

I first tried adding each of the basis functions together as function handles, e.g.:

for i = 1:N
    if i == 1
        f = @(x) c(1)*Z{1}(x);
    else
        f = @(x) f(x)+c(i)*Z{i}(x);
    end
end

Where c is a vector of the coefficients, and Z is a cell array of function handles defining each of the basis functions. However this turned out to be very slow, I assume this is because anonymous functions store the workspace that they were created in so that they can be evaluated regardless of the context in which they are called.

I next tried writing the function to an .m file explicitly. Such as:

function f = myfun(x,c)

f = c(1)*Z1+...
c(2)*Z2+...

Where Z1, Z2, are the actual basis functions (like the mess included above) hard coded in. This was faster, and further improvement was possible by using the matlabFunction command to optimize the code and precompute as much as possible. With this solution it takes ~8 sec to evaluate the function at 1e3 points.

My present solution uses the matlabFunction command to write each basis function to its own m-file, and then I write another m-file that assembles the series by calling each of the smaller basis functions. This reduced the time to evaluate 1e3 points to ~2 seconds. However, this is still way to slow for accurate integration. The domain of the function is 3 dimensional, and I'm using simp3D from the file exchange to do the integration, as its given the best speed vs. accuracy tradeoff for all of the integrators that I've tried, such as, e.g. triplequad, and mcint also from the file exchange).

SO HERE's THE QUESTION: Any ideas on how to more efficiently evaluate a large harmonic expansion like this?

p.s. The reason that the basis functions are so complicated is because the functions that I'm approximating have some interesting symmetry, so the basis functions have been symmetrized to reflect the known symmetry of the actual function.