How can I speed up this code?

I really need to speed up the code below. It is part of a stress analysis module. The idea is that the more basis polynomial terms that are taken (defined by the m, n variables in the code), the more accurate the result, i.e. the solution - which is not presented here - will converge.

When playing around with code, you can try increasing m,n for example to 4 each, 5 each, etc. You will see that computation times go up very fast.

The major bottleneck seems to be the symbolic integration.

It should also be noted that the MATLAB command 'double' in the code is very fast for m, n at around 3. However, when m,n are >5, double starts to also require significant computation. I'm not sure why.

You can propose any changes you like to this code, as my basis functions are polynomials, I do wonder whether the matlab poly functions cannot somehow speed up the calculation. Also, I do not really care if symbolic or numerical integration is used, just that it is somewhat accurate and above all, fast...

My own attempt at speeding up the computation using parfor is included at the bottom of the code, and despite roughly halving the comnputation time for large m,n, it is not nearly fast enough.

If you want the code to run very fast to test your changes/improvements quickly, then you can set m and n = 1.

Any suggestions would be highly appreciated.

%Warning: this script will clear and therefore erase your stored variables. If you do not
%desire this, delete the next two lines of code, but beware that the
%script may not function properly.
clear;
clc;

%Define geometry
L=84e-3; %length
W=24e-3; %width

%Define the elasticity tensor
C1=[303973359.614198 89929777.3924441 0 0 0 0 0 0;89929777.3924441 303973359.614198 0 0 0 0 0 0;0 0 107021791.110877 0 0 0 0 0;0 0 0 497.738940473395 80.7354743011537 -15.9810696297077 0 0;0 0 0 80.7354743011537 391.198476275344 -15.9810696297077 0 0;0 0 0 -15.9810696297077 -15.9810696297077 103.524825925731 0 0;0 0 0 0 0 0 20440000 0;0 0 0 0 0 0 0 20440000];

%Specify the degree of the approximating basis polynomials in x/y
m=3;
n=3;

%Construct basis polynomials
syms x y z
B11 = sym(zeros(m+1, 1));
B21 = sym(zeros(m+1, 1));
B12 = sym(zeros(n+1, 1));
B22 = sym(zeros(n+1, 1));

for i = 0:m
    A=factorial(m)/(factorial(i)*factorial(m-i))*(x/L)^i*(1-x/L)^(m-i);
    B11(i+1)=A;
end

for i = 0:n
    B=factorial(n)/(factorial(i)*factorial(n-i))*(y/W)^i*(1-y/W)^(n-i);
    B12(i+1)=B;
end

%Construct displacement vectors using the basis polynomials in x/y
zero_vector=zeros((m+1)*(n+1),1);
counter=0;
Vx1 = sym(zeros((m+1)*(n+1),1));
Vy1 = sym(zeros((m+1)*(n+1),1));
Vz1 = sym(zeros((m+1)*(n+1),1));
Vrx1 = sym(zeros((m+1)*(n+1),1));
Vry1 = sym(zeros((m+1)*(n+1),1));
for i = 0:m
    for j=0:n
        counter=counter+1;
        C=B11(i+1)*B12(j+1);
        Vx1(counter,1)=C; %Bernstein bases for x-displacement
        Vy1(counter,1)=C; %Bernstein bases for y-displacement
        Vz1(counter,1)=C; %Bernstein bases for z-displacement
        Vrx1(counter,1)=C; %Bernstein bases for x-rotation
        Vry1(counter,1)=C; %Bernstein bases for y-rotation
    end
end
Vx1=[Vx1;zero_vector;zero_vector;zero_vector;zero_vector]; %Displacement in x
Vy1=[zero_vector;Vy1;zero_vector;zero_vector;zero_vector]; %Displacement in y
Vz1=[zero_vector;zero_vector;Vz1;zero_vector;zero_vector]; %Displacement in z
Vrx1=[zero_vector;zero_vector;zero_vector;Vrx1;zero_vector]; %Rotation in x
Vry1=[zero_vector;zero_vector;zero_vector;zero_vector;Vry1]; %Rotation in y

%Define the strain tensor
Vxx1=diff(Vx1,x); 
Vyy1=diff(Vy1,y);
Vxy1=diff(Vx1,y);
Vyx1=diff(Vy1,x);
Vzx1=diff(Vz1,x);
Vzy1=diff(Vz1,y);
Vrxx1=diff(Vrx1,x);
Vryy1=diff(Vry1,y);
Vrxy1=diff(Vrx1,y);
Vryx1=diff(Vry1,x);

B_epsilon_1=[Vxx1.';Vyy1.';Vxy1.'+Vyx1.']; %in-plane strains
B_kappa_1=[Vrxx1.';Vryy1.';Vrxy1.'+Vryx1.']; %curvatures
B_gamma_1=[Vzy1.'-Vry1.';Vzx1.'-Vrx1.']; %shear strains
B_1=[B_epsilon_1;B_kappa_1;B_gamma_1]; %overall strain tensor

%Calculate stiffness matrix
%**This is where the longest part of the code is**
tic
K=int(int(B_1.'*C1*B_1,x,0,L),y,0,W); %Perform area integral over geometry
K=double(K); %Convert the matrix obtained to doubles
toc
%**This is where the longest part of the code ends**

%My attempt at speeding up the bottleneck stiffness matrix calculation
if matlabpool('size') == 0
    matlabpool open 4 %Assign 4 cores to task
end
tic
K1_unintegrated=triu(B_1.'*C1*B_1);
K1_unintegrated = reshape(K1_unintegrated,1,(5*(m+1)*(n+1))^2);
Batch_length=(5*(m+1)*(n+1))^2;
addendum=0;
while mod(Batch_length,4)>0
    addendum=addendum+1;
    Batch_length=Batch_length+1;
end
Batch_in=sym(zeros(4,Batch_length/4));
Batch_in(1,:)=K1_unintegrated(1:Batch_length/4);
Batch_in(2,:)=K1_unintegrated(Batch_length/4+1:2*Batch_length/4);
Batch_in(3,:)=K1_unintegrated(2*Batch_length/4+1:3*Batch_length/4);
Batch_in(4,1:Batch_length/4-addendum)=K1_unintegrated(3*Batch_length/4+1:4*Batch_length/4-addendum);
Batch_out=sym(zeros(4,Batch_length/4));
parfor i=1:4
Batch_out(i,:)=int(int(Batch_in(i,:),x,0,L),y,0,W);
end
K_fast =[Batch_out(1,:),Batch_out(2,:),Batch_out(3,:),Batch_out(4,1:Batch_length/4-addendum)];
K_fast = reshape(K_fast,5*(m+1)*(n+1),5*(m+1)*(n+1));
K_fast = K_fast+triu(K_fast,1).';
K_fast=double(K_fast);
toc