Any comment to speed up the sum calculation of a function which includes legendre polynomials?

Any comment to speed up the sum calculation of a function which includes legendre polynomials?

f(x,y)=sum(sum(legendreP(i, x)*legendreP(j, y), i=1..I), j=1..J)

ZZ=eval(f, XX, YY)

clear   
M=1000;
I=11;
J=11;
[XX,YY]=deal(linspace(-1,1,M));
    [~,n]=size(XX);
    ZZ=zeros(n);
    for k=1:n
        for l=1:n
            for i=1:I
                for j=1:J
                    ZZ(k,l) =  ZZ(k,l) +  legendreP(i, XX(k))*legendreP(j, YY(l));
                end
             end
        end 
    end

Answer 1

[Torsten]

M = 1000;
x = linspace(-1,1,M).';
y = x.';
for i = 1:11
    Lpix(:,i) = legendreP(i,x);
    Lpjy(i,:) = legendreP(i,y);
end
ZZ = zeros(numel(x),numel(y));
for i = 1:11
    for j = 1:11
        ZZ = ZZ + Lpix(:,i).*Lpjy(j,:);
    end
end

[Mehdi]

Nice!