store alternating variable for for loop

I have this code that calulates Be(cidx) for each element by calculating B four times at 4 combinations of S1 and S2. Hence in the program B is calculated 16 times, 4 times for each element. I simply want to store B at each of the 16 times for each Be(cidx). I marked B with an annotation below in the code

E = 200*10^9;
v = 0.3;
t = 0.002;
f = 4000;
angle = 0;
%% D Matrix
D = (E/(1-v^2))*[1,v,0;v,1,0;0,0,0.5-0.5*v];
%% for loop
Ke1 = zeros (16,16);
Ke2 = zeros (16,16);
Ke3 = zeros (16,16);
Ke4 = zeros (16,16);
Be1 = zeros (3,16);
Be2 = zeros (3,16);
Be3 = zeros (3,16);
Be4 = zeros (3,16);
c1 = [0,0;0.02,0;0.02,0.01;0.01,0.01;0.01,0;0.02,0.005;0.015,0.01;0.005,0.005];
c2 = [0.02,0;0.03,0;0.03,0.01;0.02,0.01;0.025,0;0.03,0.005;0.025,0.01;0.02,0.005];
c3 = [0,0;0.01,0.01;0.01,0.025;0,0.025;0.005,0.005;0.01,0.0175;0.005,0.025;0,0.0125];
c4 = [0,0.025;0.01,0.025;0.01,0.04;0,0.04;0.005,0.025;0.01,0.0325;0.005,0.04;0,0.0325];
c = {c1, c2, c3, c4};
Ke = {Ke1, Ke2, Ke3, Ke4};
Be = {Be1, Be2, Be3, Be4};
for S1 = [-1/3^0.5,1/3^0.5];
    for S2 = [-1/3^0.5,1/3^0.5];
        dN1dS1 = (((1-S2)*(S2+2*S1))/4);
        dN1dS2 = (((1-S1)*(S1+2*S2))/4);
        dN2dS1 = (((1-S2)*(2*S1-S2))/4);
        dN2dS2 = (((1+S1)*(2*S2-S1))/4);
        dN3dS1 = (((1+S2)*(2*S1+S2))/4);
        dN3dS2 = (((1+S1)*(2*S2+S1))/4);
        dN4dS1 = (((1+S2)*(2*S1-S2))/4);
        dN4dS2 = (((1-S1)*(2*S2-S1))/4);
        dN5dS1 = ((-S1*(1-S2)));
        dN5dS2 = ((-1*(1-S1^2))/2);
        dN6dS1 = ((1-S2^2)/2);
        dN6dS2 = (-S2*(1+S1));
        dN7dS1 = (-S1*(1+S2));
        dN7dS2 = ((1-S1^2)/2);
        dN8dS1 = ((-1*(1-S2^2))/2);
        dN8dS2 = (-S2*(1-S1));
        H = 1;        
        N1 = [dN1dS1,dN2dS1,dN3dS1,dN4dS1,dN5dS1,dN6dS1,dN7dS1,dN8dS1;dN1dS2,dN2dS2,dN3dS2,dN4dS2,dN5dS2,dN6dS2,dN7dS2,dN8dS2]; 
        B1 = [1,0,0,0;0,0,0,1;0,1,1,0];
        B3 = [dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1,0;dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2,0;0,dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1;0,dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2];
        
        for cidx = 1 : 4 
               J = N1*c{cidx};
               Jdet = det(J);
               B2 = [J(2,2)/Jdet,-J(1,2)/Jdet,0,0;-J(2,1)/Jdet,J(1,1)/Jdet,0,0;0,0,J(2,2)/Jdet,-J(1,2)/Jdet;0,0,-J(2,1)/Jdet,J(1,1)/Jdet];
               B = B1*B2*B3; %% alternating variable - I want to store 16 times (four times for each Be (cidx))
               Be{cidx} = Be{cidx}+ B;
               k = H*t*B'*D*B*Jdet; 
               Ke{cidx} = Ke{cidx}+k;
        end
    end
end
Ke1 = Ke{1};
Ke2 = Ke{2};
Ke3 = Ke{3};
Ke4 = Ke{4};
Be1 = Be{1};
Be2 = Be{2};
Be3 = Be{3};
Be4 = Be{4};
%% Global Stifness Matrix Assembly
KGe1 = zeros(46,46);
KGe1 ([1,2,7,8,13,14,15,16,21,22,23,24,25,26,27,28],[1,2,7,8,13,14,15,16,21,22,23,24,25,26,27,28]) = Ke1 ([1,2,7,8,3,4,5,6,9,10,11,12,13,14,15,16],[1,2,7,8,3,4,5,6,9,10,11,12,13,14,15,16]);
KGe2 = zeros(46,46);
KGe2 ([3,4,5,6,13,14,15,16,23,24,29,30,31,32,33,34],[3,4,5,6,13,14,15,16,23,24,29,30,31,32,33,34]) = Ke2 ([3,4,5,6,1,2,7,8,15,16,9,10,11,12,13,14],[3,4,5,6,1,2,7,8,15,16,9,10,11,12,13,14]);
KGe3 = zeros(46,46);
KGe3 ([1,2,7,8,17,18,19,20,27,28,35,36,37,38,39,40],[1,2,7,8,17,18,19,20,27,28,35,36,37,38,39,40]) = Ke3 ([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16],[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]);
KGe4 = zeros(46,46);
KGe4 ([9,10,11,12,17,18,19,20,37,38,41,42,43,44,45,46],[9,10,11,12,17,18,19,20,37,38,41,42,43,44,45,46]) = Ke4 ([5,6,7,8,3,4,1,2,9,10,11,12,13,14,15,16],[5,6,7,8,3,4,1,2,9,10,11,12,13,14,15,16]);
KGlobal = KGe1+KGe2+KGe3+KGe4;
%% Gaussian Elimination
KBC = zeros(40,40);
KBC ([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40],[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40]) = KGlobal ([1,2,3,4,5,6,7,8,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,45,46],[1,2,3,4,5,6,7,8,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,45,46]);
F = [0;0;(f/6)*cosd(angle);(f/6)*sind(angle);(f/6)*cosd(angle);(f/6)*sind(angle);0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;(2*f/3)*cosd(angle);(2*f/3)*sind(angle);0;0;0;0;0;0;0;0;0;0;0;0];
U = inv(KBC)*F;
%% Post Processing
%% Element 1
%% Coordinates
x1 = 0; % node (1 - local / 1 - global)
y1 = 0; % node (1 - local / 1 - global)
x7 = 0.02; % node (2 - local / 7 - global)
y7 = 0; % node (2 - local / 7 - global)
x8 = 0.02; % node (3 - local / 8 - global)
y8 = 0.01; % node (3 - local / 8 - global)
x4 = 0.01; % node (4 - local / 4 - global)
y4 = 0.01; % node (4 - local / 4 - global)
x11= 0.01; % node (5 - local / 11 - global)
y11 = 0; % node (5 - local / 11 - global)
x12 = 0.02; % node (6 - local / 12 - global)
y12 = 0.005; % node (6 - local / 12 - global)
x13 = 0.015; % node (7 - local / 13 - global)
y13 = 0.01; % node (7 - local / 13 - global)
x14 = 0.005; % node (8 - local / 14 - global)
y14 = 0.005; % node (8 - local / 14 - global)
%% Element 2
%% Coordinates
x7 = 0.02; % node (1 - local / 7 - global)
y7 = 0; % node (1 - local / 7 - global)
x2 = 0.03; % node (2 - local / 2 - global)
y2 = 0; % node (2 - local / 2 - global)
x3 = 0.03; % node (3 - local / 3 - global)
y3 = 0.01; % node (3 - local / 3 - global)
x8 = 0.02; % node (4 - local / 8 - global)
y8 = 0.01; % node (4 - local / 8 - global)
x15 = 0.025; % node (5 - local / 15 - global)
y15 = 0; % node (5 - local / 15 - global)
x16 = 0.03; % node (6 - local / 16 - global)
y16 = 0.005; % node (6 - local / 16 - global)
x17 = 0.025; % node (7 - local / 17 - global)
y17 = 0.01; % node (7 - local / 17 - global)
x12 = 0.02; % node (8 - local / 12 - global)
y12 = 0.005; % node (8 - local / 12 - global)
%% Element 3
%% Coordinates
x1 = 0; % node (1 - local / 1 - global)
y1 = 0; % node (1 - local / 1 - global)
x4 = 0.01; % node (2 - local / 4 - global)
y4 = 0.01; % node (2 - local / 4 - global)
x9 = 0.01; % node (3 - local / 9 - global)
y9 = 0.025; % node (3 - local / 9 - global)
x10 = 0; % node (4 - local / 10 - global)
y10 = 0.025; % node (4 - local / 10 - global)
x14 = 0.005; % node (5 - local / 14 - global)
y14 = 0.005; % node (5 - local / 14 - global)
x18 = 0.01; % node (6 - local / 18 - global)
y18 = 0.0175; % node (6 - local / 18 - global)
x19 = 0.005; % node (7 - local / 19 - global)
y19 = 0.025; % node (7 - local / 19 - global)
x20 = 0; % node (8 - local / 20 - global)
y20 = 0.0125; % node (8 - local / 20 - global)
%% Element 4
%% Coordinates
x10 = 0; % node (1 - local / 10 - global)
y10 = 0.025; % node (1 - local / 10 - global)
x9 = 0.01; % node (2 - local / 9 - global)
y9 = 0.025; % node (2 - local / 9 - global)
x5 = 0.01; % node (3 - local / 5 - global)
y5 = 0.04; % node (3 - local / 5 - global)
x6 = 0; % node (4 - local / 6 - global)
y6 = 0.04; % node (4 - local / 6 - global)
x19 = 0.005; % node (5 - local / 19 - global)
y19 = 0.025; % node (5 - local / 19 - global)
x21 = 0.01; % node (6 - local / 21 - global)
y21 = 0.0325; % node (6 - local / 21 - global)
x22 = 0.005; % node (7 - local / 22 - global)
y22 = 0.04; % node (7 - local / 22 - global)
x23 = 0; % node (8 - local / 23 - global)
y23 = 0.0325; % node (8 - local / 23 - global)
%% Displacement for each element
%% Element 1
u1 = U(1,1); % NODE 1
v1 = U(2,1); % Node 1
u4 = U(7,1); % NODE 4
v4 = U(8,1); % Node 4
u7 = U(9,1); % NODE 7
v7 = U(10,1); % Node 7
u8 = U(11,1); % NODE 8
v8 = U(12,1); % Node 8
u11 = U(17,1); % NODE 11
v11 = U(18,1); % Node 11
u12 = U(19,1); % NODE 12
v12 = U(20,1); % Node 12
u13 = U(21,1); % NODE 13
v13 = U(22,1); % Node 13
u14 = U(23,1); % NODE 14
v14 = U(24,1); % Node 14
%% Strain & Stress
U_e1 = [u1;v1;u7;v7;u8;v8;u4;v4;u11;v11;u12;v12;u13;v13;u14;v14];
Strain_e1 = Be1*U_e1;
Stress_e1 = D*Strain_e1;
%% Element 2
u2 = U(3,1); % NODE 2
v2 = U(4,1); % Node 2
u3 = U(5,1); % NODE 3
v3 = U(6,1); % Node 3
u7 = U(9,1); % NODE 7
v7 = U(10,1); % Node 7
u8 = U(11,1); % NODE 8
v8 = U(12,1); % Node 8
u12 = U(19,1); % NODE 12
v12 = U(20,1); % Node 12
u15 = U(25,1); % NODE 15
v15 = U(26,1); % Node 15
u16 = U(27,1); % NODE 16
v16 = U(28,1); % Node 16
u17 = U(29,1); % NODE 17
v17 = U(30,1); % Node 17
%% Strain & Stress
U_e2 = [u7;v7;u2;v2;u3;v3;u8;v8;u15;v15;u16;v16;u17;v17;u12;v12];
Strain_e2 = Be2*U_e2;
Stress_e2 = D*Strain_e2;
%% Element 3
u1 = U(1,1); % NODE 1
v1 = U(2,1); % Node 1
u4 = U(7,1); % NODE 4
v4 = U(8,1); % Node 4
u9 = U(13,1); % NODE 9
v9 = U(14,1); % Node 9
u10 = U(15,1); % NODE 10
v10 = U(16,1); % Node 10
u14 = U(23,1); % NODE 14
v14 = U(24,1); % Node 14
u18 = U(31,1); % NODE 18
v18 = U(32,1); % Node 18
u19 = U(33,1); % NODE 19
v19 = U(34,1); % Node 19
u20 = U(35,1); % NODE 20
v20 = U(36,1); % Node 20
%% Strain & Stress
U_e3 = [u1;v1;u4;v4;u9;v9;u10;v10;u14;v14;u18;v18;u19;v19;u20;v20];
Strain_e3 = Be3*U_e3;
Stress_e3 = D*Strain_e3;
%% Element 4
u5 = 0; % NODE 5
v5 = 0; % NODE 5
u6 = 0; % NODE 6
v6 = 0; % NODE 6
u9 = U(13,1); % NODE 9
v9 = U(14,1); % Node 9
u10 = U(15,1); % NODE 10
v10 = U(16,1); % Node 10
u19 = U(33,1); % NODE 19
v19 = U(34,1); % Node 19
u21 = U(37,1); % NODE 21
v21 = U(38,1); % Node 21
u22 = 0; % NODE 22
v22 = 0; % NODE 22
u23 = U(39,1); % NODE 23
v23 = U(40,1); % Node 23
%% Strain & Stress
U_e4 = [u10;v10;u9;v9;u5;v5;u6;v6;u19;v19;u21;v21;u22;v22;u23;v23];
Strain_e4 = Be4*U_e4;
Stress_e4 = D*Strain_e4;

Answer 1

[Walter Roberson]

You asked this already, and I already showed you how to do it, with exact code using cell arrays, but you deleted that question.

B = {};
for S1 = [-1/3^0.5,1/3^0.5];
    for S2 = [-1/3^0.5,1/3^0.5];
        dN1dS1 = (((1-S2)*(S2+2*S1))/4);
        dN1dS2 = (((1-S1)*(S1+2*S2))/4);
        dN2dS1 = (((1-S2)*(2*S1-S2))/4);
        dN2dS2 = (((1+S1)*(2*S2-S1))/4);
        dN3dS1 = (((1+S2)*(2*S1+S2))/4);
        dN3dS2 = (((1+S1)*(2*S2+S1))/4);
        dN4dS1 = (((1+S2)*(2*S1-S2))/4);
        dN4dS2 = (((1-S1)*(2*S2-S1))/4);
        dN5dS1 = ((-S1*(1-S2)));
        dN5dS2 = ((-1*(1-S1^2))/2);
        dN6dS1 = ((1-S2^2)/2);
        dN6dS2 = (-S2*(1+S1));
        dN7dS1 = (-S1*(1+S2));
        dN7dS2 = ((1-S1^2)/2);
        dN8dS1 = ((-1*(1-S2^2))/2);
        dN8dS2 = (-S2*(1-S1));
        H = 1;        
        N1 = [dN1dS1,dN2dS1,dN3dS1,dN4dS1,dN5dS1,dN6dS1,dN7dS1,dN8dS1;dN1dS2,dN2dS2,dN3dS2,dN4dS2,dN5dS2,dN6dS2,dN7dS2,dN8dS2]; 
        B1 = [1,0,0,0;0,0,0,1;0,1,1,0];
        B3 = [dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1,0;dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2,0;0,dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1;0,dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2];
        
        for cidx = 1 : 4 
               J = N1*c{cidx};
               Jdet = det(J);
               Bidx = 
               B2 = [J(2,2)/Jdet,-J(1,2)/Jdet,0,0;-J(2,1)/Jdet,J(1,1)/Jdet,0,0;0,0,J(2,2)/Jdet,-J(1,2)/Jdet;0,0,-J(2,1)/Jdet,J(1,1)/Jdet];
               B{end+1} = B1*B2*B3; %% alternating variable - I want to store 16 times (four times for each Be (cidx))
               Be{cidx} = Be{cidx}+ B{end};
               k = H*t*B{end}'*D*B{end}*Jdet; 
               Ke{cidx} = Ke{cidx}+k;
        end
    end
end

[Shehab Abougamrah]

Mate, I appreciate the help but I tried it before and it didn't work with the end function as it said invalid expression. Furthermore I tried this code but its missing. May you please send me the complete one based on my attached code. Cheers for all your help

[Walter Roberson]

E = 200*10^9;
v = 0.3;
t = 0.002;
f = 4000;
angle = 0;
%% D Matrix
D = (E/(1-v^2))*[1,v,0;v,1,0;0,0,0.5-0.5*v];
%% for loop
Ke1 = zeros (16,16);
Ke2 = zeros (16,16);
Ke3 = zeros (16,16);
Ke4 = zeros (16,16);
Be1 = zeros (3,16);
Be2 = zeros (3,16);
Be3 = zeros (3,16);
Be4 = zeros (3,16);
c1 = [0,0;0.02,0;0.02,0.01;0.01,0.01;0.01,0;0.02,0.005;0.015,0.01;0.005,0.005];
c2 = [0.02,0;0.03,0;0.03,0.01;0.02,0.01;0.025,0;0.03,0.005;0.025,0.01;0.02,0.005];
c3 = [0,0;0.01,0.01;0.01,0.025;0,0.025;0.005,0.005;0.01,0.0175;0.005,0.025;0,0.0125];
c4 = [0,0.025;0.01,0.025;0.01,0.04;0,0.04;0.005,0.025;0.01,0.0325;0.005,0.04;0,0.0325];
c = {c1, c2, c3, c4};
Ke = {Ke1, Ke2, Ke3, Ke4};
Be = {Be1, Be2, Be3, Be4};
B = {};
for S1 = [-1/3^0.5,1/3^0.5];
    for S2 = [-1/3^0.5,1/3^0.5];
        dN1dS1 = (((1-S2)*(S2+2*S1))/4);
        dN1dS2 = (((1-S1)*(S1+2*S2))/4);
        dN2dS1 = (((1-S2)*(2*S1-S2))/4);
        dN2dS2 = (((1+S1)*(2*S2-S1))/4);
        dN3dS1 = (((1+S2)*(2*S1+S2))/4);
        dN3dS2 = (((1+S1)*(2*S2+S1))/4);
        dN4dS1 = (((1+S2)*(2*S1-S2))/4);
        dN4dS2 = (((1-S1)*(2*S2-S1))/4);
        dN5dS1 = ((-S1*(1-S2)));
        dN5dS2 = ((-1*(1-S1^2))/2);
        dN6dS1 = ((1-S2^2)/2);
        dN6dS2 = (-S2*(1+S1));
        dN7dS1 = (-S1*(1+S2));
        dN7dS2 = ((1-S1^2)/2);
        dN8dS1 = ((-1*(1-S2^2))/2);
        dN8dS2 = (-S2*(1-S1));
        H = 1;        
        N1 = [dN1dS1,dN2dS1,dN3dS1,dN4dS1,dN5dS1,dN6dS1,dN7dS1,dN8dS1;dN1dS2,dN2dS2,dN3dS2,dN4dS2,dN5dS2,dN6dS2,dN7dS2,dN8dS2]; 
        B1 = [1,0,0,0;0,0,0,1;0,1,1,0];
        B3 = [dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1,0;dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2,0;0,dN1dS1,0,dN2dS1,0,dN3dS1,0,dN4dS1,0,dN5dS1,0,dN6dS1,0,dN7dS1,0,dN8dS1;0,dN1dS2,0,dN2dS2,0,dN3dS2,0,dN4dS2,0,dN5dS2,0,dN6dS2,0,dN7dS2,0,dN8dS2];
        
        for cidx = 1 : 4 
               J = N1*c{cidx};
               Jdet = det(J);
               B2 = [J(2,2)/Jdet,-J(1,2)/Jdet,0,0;-J(2,1)/Jdet,J(1,1)/Jdet,0,0;0,0,J(2,2)/Jdet,-J(1,2)/Jdet;0,0,-J(2,1)/Jdet,J(1,1)/Jdet];
               B{end+1} = B1*B2*B3; %% alternating variable - I want to store 16 times (four times for each Be (cidx))
               Be{cidx} = Be{cidx}+ B{end};
               k = H*t*B{end}'*D*B{end}*Jdet; 
               Ke{cidx} = Ke{cidx}+k;
        end
    end
end
Ke1 = Ke{1};
Ke2 = Ke{2};
Ke3 = Ke{3};
Ke4 = Ke{4};
Be1 = Be{1};
Be2 = Be{2};
Be3 = Be{3};
Be4 = Be{4};
%% Global Stifness Matrix Assembly
KGe1 = zeros(46,46);
KGe1 ([1,2,7,8,13,14,15,16,21,22,23,24,25,26,27,28],[1,2,7,8,13,14,15,16,21,22,23,24,25,26,27,28]) = Ke1 ([1,2,7,8,3,4,5,6,9,10,11,12,13,14,15,16],[1,2,7,8,3,4,5,6,9,10,11,12,13,14,15,16]);
KGe2 = zeros(46,46);
KGe2 ([3,4,5,6,13,14,15,16,23,24,29,30,31,32,33,34],[3,4,5,6,13,14,15,16,23,24,29,30,31,32,33,34]) = Ke2 ([3,4,5,6,1,2,7,8,15,16,9,10,11,12,13,14],[3,4,5,6,1,2,7,8,15,16,9,10,11,12,13,14]);
KGe3 = zeros(46,46);
KGe3 ([1,2,7,8,17,18,19,20,27,28,35,36,37,38,39,40],[1,2,7,8,17,18,19,20,27,28,35,36,37,38,39,40]) = Ke3 ([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16],[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]);
KGe4 = zeros(46,46);
KGe4 ([9,10,11,12,17,18,19,20,37,38,41,42,43,44,45,46],[9,10,11,12,17,18,19,20,37,38,41,42,43,44,45,46]) = Ke4 ([5,6,7,8,3,4,1,2,9,10,11,12,13,14,15,16],[5,6,7,8,3,4,1,2,9,10,11,12,13,14,15,16]);
KGlobal = KGe1+KGe2+KGe3+KGe4;
%% Gaussian Elimination
KBC = zeros(40,40);
KBC ([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40],[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40]) = KGlobal ([1,2,3,4,5,6,7,8,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,45,46],[1,2,3,4,5,6,7,8,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,45,46]);
F = [0;0;(f/6)*cosd(angle);(f/6)*sind(angle);(f/6)*cosd(angle);(f/6)*sind(angle);0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;0;(2*f/3)*cosd(angle);(2*f/3)*sind(angle);0;0;0;0;0;0;0;0;0;0;0;0];
U = inv(KBC)*F;
%% Post Processing
%% Element 1
%% Coordinates
x1 = 0; % node (1 - local / 1 - global)
y1 = 0; % node (1 - local / 1 - global)
x7 = 0.02; % node (2 - local / 7 - global)
y7 = 0; % node (2 - local / 7 - global)
x8 = 0.02; % node (3 - local / 8 - global)
y8 = 0.01; % node (3 - local / 8 - global)
x4 = 0.01; % node (4 - local / 4 - global)
y4 = 0.01; % node (4 - local / 4 - global)
x11= 0.01; % node (5 - local / 11 - global)
y11 = 0; % node (5 - local / 11 - global)
x12 = 0.02; % node (6 - local / 12 - global)
y12 = 0.005; % node (6 - local / 12 - global)
x13 = 0.015; % node (7 - local / 13 - global)
y13 = 0.01; % node (7 - local / 13 - global)
x14 = 0.005; % node (8 - local / 14 - global)
y14 = 0.005; % node (8 - local / 14 - global)
%% Element 2
%% Coordinates
x7 = 0.02; % node (1 - local / 7 - global)
y7 = 0; % node (1 - local / 7 - global)
x2 = 0.03; % node (2 - local / 2 - global)
y2 = 0; % node (2 - local / 2 - global)
x3 = 0.03; % node (3 - local / 3 - global)
y3 = 0.01; % node (3 - local / 3 - global)
x8 = 0.02; % node (4 - local / 8 - global)
y8 = 0.01; % node (4 - local / 8 - global)
x15 = 0.025; % node (5 - local / 15 - global)
y15 = 0; % node (5 - local / 15 - global)
x16 = 0.03; % node (6 - local / 16 - global)
y16 = 0.005; % node (6 - local / 16 - global)
x17 = 0.025; % node (7 - local / 17 - global)
y17 = 0.01; % node (7 - local / 17 - global)
x12 = 0.02; % node (8 - local / 12 - global)
y12 = 0.005; % node (8 - local / 12 - global)
%% Element 3
%% Coordinates
x1 = 0; % node (1 - local / 1 - global)
y1 = 0; % node (1 - local / 1 - global)
x4 = 0.01; % node (2 - local / 4 - global)
y4 = 0.01; % node (2 - local / 4 - global)
x9 = 0.01; % node (3 - local / 9 - global)
y9 = 0.025; % node (3 - local / 9 - global)
x10 = 0; % node (4 - local / 10 - global)
y10 = 0.025; % node (4 - local / 10 - global)
x14 = 0.005; % node (5 - local / 14 - global)
y14 = 0.005; % node (5 - local / 14 - global)
x18 = 0.01; % node (6 - local / 18 - global)
y18 = 0.0175; % node (6 - local / 18 - global)
x19 = 0.005; % node (7 - local / 19 - global)
y19 = 0.025; % node (7 - local / 19 - global)
x20 = 0; % node (8 - local / 20 - global)
y20 = 0.0125; % node (8 - local / 20 - global)
%% Element 4
%% Coordinates
x10 = 0; % node (1 - local / 10 - global)
y10 = 0.025; % node (1 - local / 10 - global)
x9 = 0.01; % node (2 - local / 9 - global)
y9 = 0.025; % node (2 - local / 9 - global)
x5 = 0.01; % node (3 - local / 5 - global)
y5 = 0.04; % node (3 - local / 5 - global)
x6 = 0; % node (4 - local / 6 - global)
y6 = 0.04; % node (4 - local / 6 - global)
x19 = 0.005; % node (5 - local / 19 - global)
y19 = 0.025; % node (5 - local / 19 - global)
x21 = 0.01; % node (6 - local / 21 - global)
y21 = 0.0325; % node (6 - local / 21 - global)
x22 = 0.005; % node (7 - local / 22 - global)
y22 = 0.04; % node (7 - local / 22 - global)
x23 = 0; % node (8 - local / 23 - global)
y23 = 0.0325; % node (8 - local / 23 - global)
%% Displacement for each element
%% Element 1
u1 = U(1,1); % NODE 1
v1 = U(2,1); % Node 1
u4 = U(7,1); % NODE 4
v4 = U(8,1); % Node 4
u7 = U(9,1); % NODE 7
v7 = U(10,1); % Node 7
u8 = U(11,1); % NODE 8
v8 = U(12,1); % Node 8
u11 = U(17,1); % NODE 11
v11 = U(18,1); % Node 11
u12 = U(19,1); % NODE 12
v12 = U(20,1); % Node 12
u13 = U(21,1); % NODE 13
v13 = U(22,1); % Node 13
u14 = U(23,1); % NODE 14
v14 = U(24,1); % Node 14
%% Strain & Stress
U_e1 = [u1;v1;u7;v7;u8;v8;u4;v4;u11;v11;u12;v12;u13;v13;u14;v14];
Strain_e1 = Be1*U_e1;
Stress_e1 = D*Strain_e1;
%% Element 2
u2 = U(3,1); % NODE 2
v2 = U(4,1); % Node 2
u3 = U(5,1); % NODE 3
v3 = U(6,1); % Node 3
u7 = U(9,1); % NODE 7
v7 = U(10,1); % Node 7
u8 = U(11,1); % NODE 8
v8 = U(12,1); % Node 8
u12 = U(19,1); % NODE 12
v12 = U(20,1); % Node 12
u15 = U(25,1); % NODE 15
v15 = U(26,1); % Node 15
u16 = U(27,1); % NODE 16
v16 = U(28,1); % Node 16
u17 = U(29,1); % NODE 17
v17 = U(30,1); % Node 17
%% Strain & Stress
U_e2 = [u7;v7;u2;v2;u3;v3;u8;v8;u15;v15;u16;v16;u17;v17;u12;v12];
Strain_e2 = Be2*U_e2;
Stress_e2 = D*Strain_e2;
%% Element 3
u1 = U(1,1); % NODE 1
v1 = U(2,1); % Node 1
u4 = U(7,1); % NODE 4
v4 = U(8,1); % Node 4
u9 = U(13,1); % NODE 9
v9 = U(14,1); % Node 9
u10 = U(15,1); % NODE 10
v10 = U(16,1); % Node 10
u14 = U(23,1); % NODE 14
v14 = U(24,1); % Node 14
u18 = U(31,1); % NODE 18
v18 = U(32,1); % Node 18
u19 = U(33,1); % NODE 19
v19 = U(34,1); % Node 19
u20 = U(35,1); % NODE 20
v20 = U(36,1); % Node 20
%% Strain & Stress
U_e3 = [u1;v1;u4;v4;u9;v9;u10;v10;u14;v14;u18;v18;u19;v19;u20;v20];
Strain_e3 = Be3*U_e3;
Stress_e3 = D*Strain_e3;
%% Element 4
u5 = 0; % NODE 5
v5 = 0; % NODE 5
u6 = 0; % NODE 6
v6 = 0; % NODE 6
u9 = U(13,1); % NODE 9
v9 = U(14,1); % Node 9
u10 = U(15,1); % NODE 10
v10 = U(16,1); % Node 10
u19 = U(33,1); % NODE 19
v19 = U(34,1); % Node 19
u21 = U(37,1); % NODE 21
v21 = U(38,1); % Node 21
u22 = 0; % NODE 22
v22 = 0; % NODE 22
u23 = U(39,1); % NODE 23
v23 = U(40,1); % Node 23
%% Strain & Stress
U_e4 = [u10;v10;u9;v9;u5;v5;u6;v6;u19;v19;u21;v21;u22;v22;u23;v23];
Strain_e4 = Be4*U_e4;
Stress_e4 = D*Strain_e4;

format long g
disp(Stress_e4)   %prove it got here without error

           148919229.09308
      -8.8036060333252e-05
         -800000000.000109

disp(U_e4)

      0.000422168177733788
     -0.000240830147431549
      0.000432052591607152
      0.000232093486019632
                         0
                         0
                         0
                         0
      0.000414625032753202
      3.44067134845473e-06
      0.000125042577101716
      0.000128037790548056
                         0
                         0
      0.000124721445024561
     -0.000133518117814768

[Shehab Abougamrah]

Cheers man it works, appreciate all the help