Calculate strain between elements using system of linear equation

I have 10 elements with their serial number and x,y,z coordinates in their initial and final position : data0 and data1 respectively.

I intend to calculate the strain as follows:

1) for each element in data0, I will search for 3 elements closest to it and for each of the four close elements (i, j, k, l), the displacements ui, uj, uk, ul in the x direction; vi, vj, vk, vl in the y direction; and wi, wj, wk, wl in the z direction are all calculated

2) simultaneous equations are created for each 4 element set to obtain a0, a1, a2; b0 ,b1 ,b2; c0 ,c1 ,c2

-[1, xi,yi,zi;1,xj,yj,zj;1,xk,yk,zk;1,xl,yl,zl] * [a0,a1,a2,a3] = [ui,uj,uk,ul]

-[1, xi,yi,zi;1,xj,yj,zj;1,xk,yk,zk;1,xl,yl,zl] * [b0,b1,b2,b3] = [vi,vj,vk,vl]

-[1, xi,yi,zi;1,xj,yj,zj;1,xk,yk,zk;1,xl,yl,zl] * [c0,c1,c2,c3] = [wi,wj,wk,wl]

3) finally for each element the following strain is calculated:

-Ex = du/dx = a1

-Ey = dv/dy = b2

-Ez = dw/dz = c3

-Exy = (b1 +a2)/2

-Exz = (a3 + c1)/2

-Eyz = (b3 + c2)/2

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------I tried creating the following code but got my simulataneous equation loop is not well structured

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------

clear;
data0 = [1, 46.15, 29.55, 4.36; 2, 56.19, 34.22, 5.33;3, 39.01, 33.28, 7.27; 4, 29.51, 37.41, 5.89;5, 65.51, 37.97, 5.04;6, 50.78 ,47.22, 3.25; 7, 72.69, 47.58 ,4.34; 8, 30.07, 48.18, 5.62; 9, 41.34, 50.11, 3.91;10, 61.13, 52.16, 2.22]                                                                                                                                                                                                                                  
%Extract element number and store element_number
Element_No1=data0(:,1);
X1=data0(:,2);
Y1=data0(:,3);
Z1=data0(:,4);
EUC_DIST= data0(:,5);

data1=[1,46.83,30,4.69;2,56.84,35.49,6.51;3,39.39,33.56,7.11;4,29.74,38.16,5.87;5,64.5,38.55,5.61;6,50.74,47.57,3.22;7,71.88,47.8,4.47;8,29.72,48.2,5.59;9,41.13,50.23,3.88;10,60.61,52.4,2.29]
Element_No2=data1(:,1);
X2=data1(:,2);
Y2=data1(:,3);
Z2=data1(:,4);
%Dx= data1(:,5);
%Dy= data1(:,6);
%Dz= data1(:,7);

%Check matrix size           
[m,n]=size(data0);
[p,o]=size(data1);

%count the number of data points
Q0 = numel(data0(:,1));
Q1 = numel(data1(:,1));
q1 = [1:Q1]';

%create matrix of distance difference with zeros
DIST_DIFF = zeros(Q0,Q1);

%Calculate the eucledean distsnce difference between each element and others
for j=1:m
    for i=1:m
                DIST_DIFF(j,i) = abs((data0(j,5)-data0(i,5)));
                DIST_DIFF(DIST_DIFF == 0) = NaN;
    end 
end

%add the element number to the calculated distance difference 
DIST_DIFF = horzcat(DIST_DIFF,q1);


%Get the closest 3 element counterpart and store them in EUCL_DIST
[r,s]=size(DIST_DIFF);
A = [DIST_DIFF(:,1),DIST_DIFF(:,(Q0+1))];
SA = size(A);
NA = zeros(SA);
EUCL_DIST = [];

for i= 1:s-1
    %generate matrix for each column of DIST_DIFF and Attach element number to each generated column
    NA = [DIST_DIFF(:,i),DIST_DIFF(:,(Q0+1))];
    %sort each generated column 
    S_NA= sortrows(NA);
    %take out the 3 closest element
    T= [S_NA(1:3,1:2)]; %3 here indicate number of elements to display 
    EUCL_DIST = [EUCL_DIST,T];
end


%Attach other element properties (x,y,z) to matched particles ID

%Isolate the 3 closest particle number and leave out the calculated distance
y0= EUCL_DIST(:,2:2:end);

[r,s]=size(y0);
TOG = [];

for i= 1:s
    %generate matrix for each column of DIST_DIFF and Attach element number to each generated column take out the 3 closest element
    q=ismember(data0(:,1), y0(:,i), 'rows');
    z=find(q);
    Close_Elements =data0(z,1:4);
    TOG = [TOG,Close_Elements];
end


ALL_VAR = zeros(Q0,21);

[p,t]=size(TOG);
[g,h]=size(data0);
[q,r]=size(data1);

for i= 1:t
   for ii= 1:h
       for iii= 1:r
        A1 =[1, data0(1,2), data0(1,3), data0(1,4);1,  TOG(1,2), TOG(1,3),  TOG(1,4);1, TOG(2,2), TOG(2,3),  TOG(2,4);1,      TOG(3,2), TOG(3,3),  TOG(3,4)]; 
        A2 =[1, data0(1,2), data0(1,3), data0(1,4);1,  TOG(1,2), TOG(1,3),  TOG(1,4);1, TOG(2,2), TOG(2,3),  TOG(2,4);1, TOG(3,2), TOG(3,3),  TOG(3,4)]; 
        A3 =[1, data0(1,2), data0(1,3), data0(1,4);1,  TOG(1,2), TOG(1,3),  TOG(1,4);1, TOG(2,2), TOG(2,3),  TOG(2,4);1, TOG(3,2), TOG(3,3),  TOG(3,4)]; 


        B1 = [(data0(1,2)-data1(1,2));(data0(2,2)-data1(2,2));(data0(3,2)-data1(3,2));(data0(4,2)-data1(4,2))];
        B2 = [(data0(1,3)-data1(1,3));(data0(2,3)-data1(2,3));(data0(3,3)-data1(3,3));(data0(4,3)-data1(4,3))];
        B3 = [(data0(1,4)-data1(1,4));(data0(2,4)-data1(2,4));(data0(3,4)-data1(3,4));(data0(4,4)-data1(4,4))];

        a = linsolve(A1,B1);
        b = linsolve(A2,B2);
        c = linsolve(A3,B3);

        Const = [a,b,c];
       end
   end
end
%%%%%Difference between initial and final particle position
Du = X2 - X1;
Dv= y2 - y1;
Dw= z2 - z1;
%%%Difference between each close particle and the particle under consideration
Dx= x2nd -x1;
Dy= x3rd-x1;
Dz= x4th-x1;

%Calculate micro strain - update du ,dv ,dw ,dx ,dy ,dz to ALL_VARIABLES
a1 = du/dx ;
b2 = dv/dy;
c3 = dw/dz;

Ex = a1;
Ey = b2;
Ez = c3;