1 function lambda = c0P2(h)
  2 %% Mesh
  3 [node,elem] = squarequadmesh([0,1,0,1],h);
  4 elem = elem(:,[1,4,3,2]);
  5   showmesh(node,elem);
  6   findnode(node);
  7   findquadelem(node,elem);
  8 %% Construct Data Structure
  9 [elem2dof,edge,inDof] = c0dofP2(elem);
 10 elem2dof=double(elem2dof);
 11 N = size(node,1);  NT = size(elem,1);  
 12 Ndof = N+NT+size(edge,1);
 13 A=sparse(Ndof,Ndof);
 14 B=sparse(Ndof,Ndof);
 15 %% Assemble stiffness matrix
 16 % Since Dphi_i*Dphi_j is quadratic, 
 17 % numerical quadrature rule is used here
 18 option.quadorder = 4;   % default order
 19 [pts,weight] = quadquadpts(option.quadorder);
 20 pts=pts*2-1;
 21 x=pts(:,1);y=pts(:,2);
 22 h1=h/2;h2=h1;
 23 area=h2*h1;
 24 %% Assemble Matrix 
 25 for i=1:9 
 26     for j=i:9
 27 DuDv=area*(Dxphi(x,y,h1,i).*Dxphi(x,y,h1,j)...
 28           +Dyphi(x,y,h2,i).*Dyphi(x,y,h2,j))‘*weight; 
 29 uv=area*(phi(x,y,i).*phi(x,y,j))‘*weight; 
 30       if i==j
 31 A = A + sparse(elem2dof(:,i),elem2dof(:,j),DuDv,Ndof,Ndof);
 32 B = B + sparse(elem2dof(:,i),elem2dof(:,j),uv,Ndof,Ndof);
 33       else
 34 A = A + sparse([elem2dof(:,i);elem2dof(:,j)],...
 35                [elem2dof(:,j);elem2dof(:,i)],DuDv,Ndof,Ndof);
 36 B = B + sparse([elem2dof(:,i);elem2dof(:,j)],...
 37                [elem2dof(:,j);elem2dof(:,i)],uv,Ndof,Ndof);
 38       end
 39     end
 40 end
 41 %% Solve Ax = lambda Bx, and its first solution is 2*pi^2
 42 lambda=eigs(A(inDof,inDof),B(inDof,inDof),1,‘sm‘); 
 43 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 44 % subfunction Dxphi
 45 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 46     function s = Dxphi(xi,eta,L1,i) % gradient of basis phi for x
 47     switch i
 48         case 1
 49             s = (eta.*(2.*xi - 1).*(eta - 1))/(4*L1);           
 50         case 2
 51             s = (eta.*(2.*xi + 1).*(eta - 1))/(4*L1);           
 52         case 3
 53             s = (eta.*(2.*xi + 1).*(eta + 1))/(4*L1);          
 54         case 4
 55             s = (eta.*(2.*xi - 1).*(eta + 1))/(4*L1);
 56         case 5
 57             s = -(eta.*xi.*(eta - 1))/L1;
 58         case 6
 59             s = -((eta.^2 - 1).*(2*xi + 1))/(2*L1);
 60         case 7
 61             s = -(eta.*xi.*(eta + 1))/L1;
 62         case 8
 63             s = -((eta.^2 - 1).*(2.*xi - 1))/(2*L1); 
 64         case 9
 65             s = (2*xi.*(eta.^2 - 1))/L1;
 66     end
 67     end
 68 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 69 % subfunction Dxphi
 70 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 71     function s = Dyphi(xi,eta,L2,i) % gradient of basis phi for y
 72      switch i
 73      case 1 
 74          s = (xi.*(2.*eta - 1).*(xi - 1))/(4*L2);
 75      case 2 
 76          s = (xi.*(2.*eta - 1).*(xi + 1))/(4*L2);
 77      case 3
 78          s = (xi.*(2.*eta + 1).*(xi + 1))/(4*L2);
 79      case 4
 80          s = (xi.*(2.*eta + 1).*(xi - 1))/(4*L2);
 81     case 5 
 82         s = -((2.*eta - 1).*(xi.^2 - 1))/(2*L2);
 83     case 6
 84         s = -(eta.*xi.*(xi + 1))/L2;
 85     case 7 
 86         s =-((2*eta + 1).*(xi.^2 - 1))/(2*L2);
 87     case 8 
 88         s = -(eta.*xi.*(xi - 1))/L2;
 89     case 9
 90         s = (2*eta.*(xi.^2 - 1))/L2;
 91      end
 92     end
 93 
 94 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 95 % subfunction phi
 96 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 97     function s = phi(xi,eta,i) % gradient of basis phi
 98     switch i
 99         case 1
100             s = (eta.*xi.*(eta - 1).*(xi - 1))/4;           
101         case 2
102             s = (eta.*xi.*(eta - 1).*(xi + 1))/4;            
103         case 3
104             s = (eta.*xi.*(eta + 1).*(xi + 1))/4;             
105         case 4
106             s = (eta.*xi.*(eta + 1).*(xi - 1))/4; 
107         case 5
108             s = -(eta.*(xi.^2 - 1).*(eta - 1))/2;
109         case 6
110             s = -(xi.*(eta.^2 - 1).*(xi + 1))/2;
111         case 7
112             s = -(eta.*(xi.^2 - 1).*(eta + 1))/2;
113         case 8
114             s = -(xi.*(eta.^2 - 1).*(xi - 1))/2;
115         case 9
116             s = (eta.^2 - 1).*(xi.^2 - 1);
117     end
118     end
119 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
120 end

 1 function [elem2dof,edge,inDof] = c0dofP2(elem)
 2 totalEdge=sort([elem(:,[1,2]);elem(:,[2,3]);elem(:,[3,4]);elem(:,[4,1])],2);
 3 [edge, i2, j] = myunique(totalEdge);
 4 N = max(elem(:)); 
 5 NT = size(elem,1);
 6 NE = size(edge,1);
 7 elem2edge = reshape(j,NT,4);
 8 elem2dof = uint32([elem N+elem2edge (N+NE+1:N+NE+NT)‘]);
 9 i1(j(4*NT:-1:1)) = 4*NT:-1:1; 
10 i1 = i1‘;
11 bdEdgeIdx = (i1 == i2);
12 isBdDof = false(N+NE+NT,1);
13 isBdDof(edge(bdEdgeIdx,:)) = true;   % nodal 
14 idx = find(bdEdgeIdx);
15 isBdDof(N+idx) = true;
16 inDof = find(~isBdDof);
17 end