To build assembly matrix from elemental matrix

here is the code below. i tried to assemble elemental matrix using a for loop but the result ends with a matrix having a determinant 0. i need to use iteration to solve the resulting equation but the convergence of the matrix doesn't allow. i dont really know where i make the mistake whether if it is on the code or somewhere else. can i get any default assemblying code or can anyone help to assess my code and idetify where i made a mistake. the analysis is based on finite element method.

%% global matrix

clear all
clc
G=[112 28 29;
    153 107 113;
%     150 110 108;
%     174 142 149;
    138 89 91;
%     172 170 165;
    63 161 168;
%     157 122 158;
%     175 115 154;
    66 67 33]


%% stiffness matrix

syms z e
mu=0.28
density=0.05
for i=1:3
    if i==1
        s(1)=1-z-e;
    elseif i==2
        s(2)=z;
    else
        s(3)=e;
    end
end
s(1)
s(2)
s(3)

%% jacobian matrix 

jj=[diff(s(1),z) diff(s(2),z) diff(s(3),z);
    diff(s(1),e) diff(s(2),e) diff(s(3),e)]
nodeco=[32.92 17.76;30.76 20.4;29.2 16.77]
jacobian=jj*nodeco
invJ=inv(jacobian)
D=det(jacobian)

%% stiffness matrix of velocity 

for i=1:3
    
for j=1:3
ku(i,j)=int(int((invJ(2,1)*diff(s(i),z)+invJ(2,2)*diff(s(j),e))*(invJ(2,1)*diff(s(i),z)+invJ(2,2)*diff(s(j),e)),e,0,1-z),z,0,1)*D*mu;
end
end
ku

%% stiffness for pressure

for i=1:3
    
for j=1:3
kp(i,j)=int(int(s(i)*(invJ(1,1)*diff(s(j),z)+invJ(1,2)*diff(s(j),e)+invJ(2,1)*diff(s(j),z)+invJ(2,2)*diff(s(j),e)),e,0,1-z),z,0,1)*D/density;
end
end
kp







%%assembled velocity coefficient 


sum=zeros(168);

for i=1:5
   K=zeros(168);
    for j=1:3
        
%       [M(j,1),M(j,2), M(j,3)]==[K(G(i,j),G(i,1)) K(G(i,j),G(i,2)) K(G(i,j),G(i,3))]
      K(G(i,j),G(i,1))=ku(j,1);
      K(G(i,j),G(i,2))=ku(j,2);
      K(G(i,j),G(i,3))=ku(j,3);
      
    end
   
sum=sum+K;
end
kuu=sum
 kuuu=kuu([28 29 33 63 66 67 89 91 107 112 113 138 153 161 168],[28 29 33 63 66 67 89 91 107 112 113 138 153 161 168])

%%assembled pressure coefficient 

sum=zeros(168);
for i=1:5
   K=zeros(168);
    for j=1:3
        
%       [M(j,1),M(j,2), M(j,3)]==[K(G(i,j),G(i,1)) K(G(i,j),G(i,2)) K(G(i,j),G(i,3))]
      K(G(i,j),G(i,1))=kp(j,1);
      K(G(i,j),G(i,2))=kp(j,2);
      K(G(i,j),G(i,3))=kp(j,3);
      
    end
   
sum=sum+K;
end
kpp=sum
kppp=kpp([28 29 33 63 66 67 89 91 107 112 113 138 153 161 168],[28 29 33 63 66 67 89 91 107 112 113 138 153 161 168])

%% source term

syms B Ui double
B=20;


conductivity=100;
for i=1:3
f(i,1)=int(int((s(i)*(B^2)*Ui),e,0,1-z),z,0,1)*(conductivity*D)/density;

end
f

%% assemble of source term 

sum=zeros([168,1]);
for i=1:5
    
   K=sym(zeros([168,1]))  ;
   for j=1:3
     K(G(i,j),1)=f(j,1);
    end
   
sum=sum+K;


end
kff=sum;
kfff=kff([28 29 33 63 66 67 89 91 107 112 113 138 153 161 168],1)
% double(kfff)

[Walter Roberson]

syms B Ui double

Caution: that would create a symbolic variable named double . There is no double modifier for symbolic variables. There are modifiers such as real and positive

[Walter Roberson]

sum=zeros([168,1]);

It is recommended to avoid using sum as the name of a variable, as it conflicts with the MATLAB function sum() . It is very common for people to use sum as a variable name and then a few lines further down try to use sum() as a function.