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

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

% ====================================================%
% function to create Q9 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/2/m:L,0:H/2/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:2:2*m
      elem_i = [[j j+2] [j+2 j]+(2*m+1)*2 j+1 ...
                (j+2)+(2*m+1) (j+1)+(2*m+1)*2 j+(2*m+1) (j+1)+(2*m+1)]...
                + (2*m+1)*2*(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.1)
         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