Create Q4 Elements for 2d Rectangular Domain (Octave/MATLAB)

function [node_xy,elem_node] = ElemQ4(L,H,m,n,varargin)

% ====================================================%
% function to create Q4 elements (2D)                 %
% INPUT:                                              %
% > length (L) and height (H),                        %
% > division over x (m) and division over y (n)       %
% > option for plotting (varargin)                    %
% OUTPUT:                                             %
% > coordinate of nodes [node_xy]                     %
% > element connectivity [elem_node]                  %
% ====================================================%

% Nodes coordinates [x y]
[x,y] = meshgrid(0:L/m:L,0:H/n:H);
node_xy = [x(:) y(:)];
node_xy = sortrows(node_xy,2);

% Connectivity matrix [node1 node2 node3 node4]
elem_node = [];
for i = 1:n
   for j = 1:m
      elem_i = [[j j+1] [j+1 j]+(m+1)] + (m+1)*(i-1);
      elem_node = cat(1,elem_node,elem_i);
   end
end


%===================================
% Plot (optional)
%===================================
if nargin > 4
   figure
   plot(node_xy(:,1),node_xy(:,2),'ro')
   hold on 
   for j = 1:size(elem_node,1)
      for k = 1:size(elem_node,2)-1
         pause(0.2)
         plot(node_xy(elem_node(j,k:k+1),1),...
         node_xy(elem_node(j,k:k+1),2),'b:','linewidth',3)
      end
   end 
   hold off
end

end