nested loop optimisation

I'm really struggling to optimize a calculus code on MATLAB.

It's a heavy calculation required to obtain the material properties of a non-linear material.

This calculation requires more than 240 millions steps. It's rather simple in itself as it consist in a huge number of sum. The only problem is that the numbers are stored in various arrays and lists which is a bit confusing.

Here is the original code :

Tensor=zeros(3,3,3,3);
for m=1:3
    for n=1:3
        for o=1:3
            for p=1:3
                for x=1:16 
                    for y=1:16
                        for z=1:16 
                            for i=1:3
                                for j=1:3
                                    for k=1:3
                                        for l=1:3
                                            for r=1:3
                                                for s=1:3
                                                    sum=sum+(1/(8*(pi()^2))*P{x,y,z}(i,m)*P{x,y,z}(j,n)*P{x,y,z}(k,o)*P{x,y,z}(l,p)*(TensorC(i,j,k,l)-TensorC0(i,j,r,s))*TensorA(r,s,k,l)*sin(omega(x))*p_omega(x)*p_phi(y)*p_beta(z);
                                                end
                                            end
                                        end                                                    
                                    end
                                end
                            end
                        end
                    end
                end
                Tensor(m,n,o,p)=sum;
            end
        end
    end
end

P{x,y,z}(i,m) is a change of basis formula (same for the others) : i and m determine the type of formula and the result is calculated with the x, y and z variables.

All the other numbers in the summation are real numbers picked up in arrays and tensors.

I tried to extract as much variables from the last for loop in order to calculate them as soon as possible an reduce the number of operations:

Tensor=zeros(3,3,3,3);
CO1=1/(8*(pi()^2));
for m=1:3
    for n=1:3
        for o=1:3
            for p=1:3
            sum=C0_tensor(m,n,o,p);
                for x=1:16
                    CO7=sin(omega(x));
                    CO8=p_omega(x);
                    for y=1:16
                        CO9=p_phi(y);
                        for z=1:16
                            CO10=p_beta(z);
                            for i=1:3
                                CO2=P{x,y,z}(i,m);
                                for j=1:3
                                    CO3=P{x,y,z}(j,n);
                                    for k=1:3
                                        CO4=P{x,y,z}(k,o);
                                        for l=1:3
                                            CO5=P{x,y,z}(l,p);
                                            CO6=TensorC(i,j,k,l);
                                            for r=1:3
                                                for s=1:3
                                                    CO11=TensorC0(i,j,r,s);
                                                    CO12=TensorA(r,s,k,l);
                                                    sum=sum+CO1*CO2*CO3*CO4*CO5*(CO6-CO11)*CO12*CO7*CO8*CO9*CO10;
                                                end
                                            end
                                        end                                                    
                                    end
                                end
                            end
                        end
                    end
                end
                Tensor(m,n,o,p)=sum;
            end
        end
    end
end

Still, the calculation is much too long.

I don't see any way of parallelizing or vectorizing the calculation ...

The operation of retrieving one particular value from one array or one matrix seem to be very slow...

Do you think I should build a huge tensor containing all the values instead of using multiples one?