Specifying an element and gaussian point for a script

I was wondering how exactly I would specify an individual element and gaussian point for the following:

function res = FEMbmatrix(L,n) %inputs n[2,2]=4 i.e. number of plates. L is length (2,2] squr

nx = n(1) + 1; % number of nodes = number of elements +1 (x) ;num. of nodes given num. element input established
ny = n(2) + 1; % number of nodes = number of elements +1 ? ; n will be a 1*2 matrix; n = [x,y]
nNodes = nx * ny; % total number of nodes
nE = n(1) * n(2); % total number of elements; e.g n[2,2]=4 elements
dx = L(1)/n(1); % length of x/number of elements in x
dy = L(2)/n(2); % height/number of element in y direction

% Gaussian points
Zeta = [-(1/3)^0.5; -(1/3)^0.5; (1/3)^0.5; (1/3)^0.5]; % local zeta; bottom left quadrant to bottom right quadrant clockwise.
Eta = [-(1/3)^0.5; (1/3)^0.5; (1/3)^0.5; -(1/3)^0.5]; % local eta;

eID = 1; %loop over all elemets
for j = 1:n(2)
    for i = 1:n(1)
        %node ID
        e(eID).n(1) = (j-1)*nx+i;
        e(eID).n(2) = j*nx+i;
        e(eID).n(3) = j*nx+i+1;
        e(eID).n(4) = (j-1)*nx+i+1;
        %global node position
        e(eID).x(1,1) = (i-1)*dx;
        e(eID).y(1,1) = (j-1)*dy;
        e(eID).x(2,1) = (i-1)*dx;
        e(eID).y(2,1) = j*dy;
        e(eID).x(3,1) = i*dx;
        e(eID).y(3,1) = j*dy;
        e(eID).x(4,1) = i*dx;
        e(eID).y(4,1) = (j-1)*dy;
        eID=eID+1;
    end
end


for eID = 1:nE
   
    for gp=1:4                 %loop over all gaussian points
        zeta = Zeta(gp,1);
        eta = Eta(gp,1);
        dNdzeta(1) = -0.25*(1-eta);
        dNdzeta(2) = -0.25*(1+eta);
        dNdzeta(3) = 0.25*(1+eta);
        dNdzeta(4) = 0.25*(1-eta);
        dNdeta(1) = -0.25*(1-zeta);
        dNdeta(2) = 0.25*(1-zeta);
        dNdeta(3) = 0.25*(1+zeta);
        dNdeta(4) = -0.25*(1+zeta);
        dNdlocal = [dNdzeta; dNdeta];
        %Jacobian
        J(1,1) = dNdzeta*e(eID).x;
        J(1,2) = dNdzeta*e(eID).y;
        J(2,1) = dNdeta*e(eID).x;
        J(2,2) = dNdeta*e(eID).y;
        dNdglobal = J\dNdlocal; %compacted matrix 
        %define 3*8 b matrix
        B(1,1) = dNdglobal(1,1);
        B(1,2) = 0;
        B(1,3) = dNdglobal(1,2);
        B(1,4) = 0;
        B(1,5) = dNdglobal(1,3);
        B(1,6) = 0;
        B(1,7) = dNdglobal(1,3);
        B(1,8) = 0;
        B(2,1) = 0;
        B(2,2) = dNdglobal(2,1);
        B(2,3) = 0;
        B(2,4) = dNdglobal(2,2);
        B(2,5) = 0;
        B(2,6) = dNdglobal(2,3);
        B(2,7) = 0;
        B(3,1) = dNdglobal(2,1);
        B(3,2) = dNdglobal(1,1);
        B(3,3) = dNdglobal(2,2);
        B(3,4) = dNdglobal(1,2);
        B(3,5) = dNdglobal(2,3);
        B(3,6) = dNdglobal(1,3);
        B(3,7) = dNdglobal(2,4);
        B(3,8) = dNdglobal(1,4);
        B    
    
    end
   
end

I have written this script for the B matrix for each point. Ideally, I would want to store each iteration of the gaussian point for each element somewhere.

Thanks!