Structured mesh in MATLAB PDE toolbox

I found a way to create an structured quadratic mesh. They way I do it is as follows. For example if I have a rectangle. I divide it to many squares and then mesh each square to two triangular mesh.

if true
numberOfPDE = 3;
model = createpde(numberOfPDE);
global dl
L=2;
h=1;
dl=.05 ;
dh=dl ;
[xq,yq]=meshgrid( [0:dl:L], [-h/2:dh:h/2]);
n=size(xq,1);
m=size(xq,2);
k=1;
R=[];
for i=1:n-1
for j=1:m-1
 R(:,k) =[3,4,xq(i,j),xq(min(i+1,n),j),xq(min(i+1,n),min(m,j+1)),xq(i,min(m,j+1)),yq(i,j),yq(min(i+1,n),j),yq(min(i+1,n),min(m,j+1)),yq(i,min(m,j+1))];    
  k=k+1;
end
end
end
gm = [R];
for i=1:size(R,2)
  str{i}  =  ( sprintf('R%d',i) );
end
sf=[];
for i=1:size(R,2)
  sf=char([str{i},'+',sf]);
end
sf=sf(1:end-1);
delete(model.Geometry)
model.Geometry = []
ns = char(str);
ns = ns';
g=  decsg(gm,sf,ns); 
geometryFromEdges(model,g);
figure
pdegplot(model,'EdgeLabels','on');
axis equal
title 'Geometry with Edge Labels';

This code creates the geometry, The following is the meshing part.

if true
order='quadratic';
% order='linear';
mesh=generateMesh(model,'GeometricOrder', order    ,'Hmax' ,.05  , 'MesherVersion','R2013a','Hgrad',1.01,'JiggleIter',100  );
figure(200)
pdemesh(model,'NodeLabels','off')
[p,e,t] = meshToPet(mesh);
axis equal
end

Then PDE is solved as follows.

if true
a = [0 0 0 0 0 0 0 0 0]';
f = [0 0 0]’;
applyBoundaryCondition(model,'mixed',...
  'Edge',  761:780   ,...
  'u',  0,'EquationIndex',  3  ,'Vectorized','on'  );
applyBoundaryCondition(model,'dirichlet',...
  'Edge',[1:40,1601:1609],...
  'u',  @myConstPsi2Rev,'Vectorized','on'  );
applyBoundaryCondition(model,'dirichlet',...
  'Edge',1610:1611,...
  'u',  @myConstPsi3Rev,'Vectorized','on'  );
applyBoundaryCondition(model,'dirichlet',...
  'Edge',[1612:1620,1561:1600],...
  'u',  @myConstPsi4Rev,'Vectorized','on'  );
C=@ccoefPsiDevi;
specifyCoefficients(model,'m',0,'d',0,'c',C,'a',a,'f',f);
end

if true
  format long
model.SolverOptions.MinStep=0;
model.SolverOptions.ReportStatistics = 'on';
model.SolverOptions.MaxIterations=5000;
model.SolverOptions.ResidualTolerance=1.6819e-06
%%solve PDE
tic
result = solvepde(model);
toc
end

The functions used as C coefficients and Boundary conditions are as follows.

if true
function cmatrix  = ccoefPsiDevi(region,state)
%   global xe ye  %epsilon %error
n1 = 36;
nr = numel(region.x) ;
cmatrix = zeros(n1,nr);
  e= (.1);
  normrev= (sqrt(e^2 +(state.uy(1,:)+state.ux(2,:)).^2+2*state.uy(2,:).^2+2*state.ux(1,:).^2     ) );
fac=1 ;
cmatrix(1,:) =   fac*2./normrev;
cmatrix(4,:) =   fac*1./normrev;
cmatrix(14,:) =   fac*1./normrev;
cmatrix(26,:) =  -region.y;
cmatrix(27,:) =   region.y;
cmatrix(7,:) =   fac*1./normrev;
cmatrix(17,:) =  fac*1./normrev;
cmatrix(20,:) =  fac*2./normrev;
cmatrix(30,:) =  region.x;
cmatrix(31,:) = -region.x;
cmatrix(10,:) =  -region.y;
cmatrix(11,:) =   region.y;
cmatrix(22,:) =  region.x;
cmatrix(23,:) = -region.x;            
end
end

if true
function bcMatrix = myConstPsi2Rev(region,state)
bcMatrix = [    0*region.y ;  
            20*dl*ones(size(region.y));    
            0*region.y] ;
end
end

if true
function bcMatrix = myConstPsi3Rev(region,state)
global psixb psiyb  
bcMatrix = [ 0*region.y;           
-20*region.y;
0*region.y;]
end
end

if true
function bcMatrix = myConstPsi4Rev(region,state)
bcMatrix = [ 0*region.y;            
-20*dl*ones(size(region.y));
           0*region.y ];
end
end

But the problem is because the geometry is consisting of many squares instead of one rectangle. My code becomes very slow. How can I change the geometry after meshing to the rectangle with out having any issues. Or is there any other way to do this. I can not use legacy work flow because the elements are only linear in legacy, and I need quadratic elements. My results are much smother with this structured mesh, rather than the Matlab default unstructured mesh.So I really need this structured mesh, but because of the geometry the code is super slow. Please help me.

[Kaveh Gharibi]

The soloution to this problem is here:

https://www.mathworks.com/matlabcentral/answers/346410-is-there-a-way-to-import-2d-meshes-nodes-and-triangles-into-pde-model-objects