function lambda = c0P2(h)

 %% Mesh

 [node,elem] = squarequadmesh([,,,],h);

 elem = elem(:,[,,,]);

   showmesh(node,elem);

   findnode(node);

   findquadelem(node,elem);

 %% Construct Data Structure

 [elem2dof,edge,inDof] = c0dofP2(elem);

 elem2dof=double(elem2dof);

 N = size(node,);  NT = size(elem,);

 Ndof = N+NT+size(edge,);

 A=sparse(Ndof,Ndof);

 B=sparse(Ndof,Ndof);

 %% Assemble stiffness matrix

 % Since Dphi_i*Dphi_j is quadratic,

 % numerical quadrature rule is used here

 option.quadorder = ;   % default order

 [pts,weight] = quadquadpts(option.quadorder);

 pts=pts*-;

 x=pts(:,);y=pts(:,);

 h1=h/;h2=h1;

 area=h2*h1;

 %% Assemble Matrix

 for i=:

     for j=i:

 DuDv=area*(Dxphi(x,y,h1,i).*Dxphi(x,y,h1,j)...

           +Dyphi(x,y,h2,i).*Dyphi(x,y,h2,j))'*weight;

 uv=area*(phi(x,y,i).*phi(x,y,j))'*weight;

       if i==j

 A = A + sparse(elem2dof(:,i),elem2dof(:,j),DuDv,Ndof,Ndof);

 B = B + sparse(elem2dof(:,i),elem2dof(:,j),uv,Ndof,Ndof);

       else

 A = A + sparse([elem2dof(:,i);elem2dof(:,j)],...

                [elem2dof(:,j);elem2dof(:,i)],DuDv,Ndof,Ndof);

 B = B + sparse([elem2dof(:,i);elem2dof(:,j)],...

                [elem2dof(:,j);elem2dof(:,i)],uv,Ndof,Ndof);

       end

     end

 end

 %% Solve Ax = lambda Bx, and its first solution is *pi^

 lambda=eigs(A(inDof,inDof),B(inDof,inDof),,'sm');

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 % subfunction Dxphi

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     function s = Dxphi(xi,eta,L1,i) % gradient of basis phi for x

     switch i

         case

             s = (eta.*(.*xi - ).*(eta - ))/(*L1);

         case

             s = (eta.*(.*xi + ).*(eta - ))/(*L1);

         case

             s = (eta.*(.*xi + ).*(eta + ))/(*L1);

         case

             s = (eta.*(.*xi - ).*(eta + ))/(*L1);

         case

             s = -(eta.*xi.*(eta - ))/L1;

         case

             s = -((eta.^ - ).*(*xi + ))/(*L1);

         case

             s = -(eta.*xi.*(eta + ))/L1;

         case

             s = -((eta.^ - ).*(.*xi - ))/(*L1);

         case

             s = (*xi.*(eta.^ - ))/L1;

     end

     end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 % subfunction Dxphi

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     function s = Dyphi(xi,eta,L2,i) % gradient of basis phi for y

      switch i

      case

          s = (xi.*(.*eta - ).*(xi - ))/(*L2);

      case

          s = (xi.*(.*eta - ).*(xi + ))/(*L2);

      case

          s = (xi.*(.*eta + ).*(xi + ))/(*L2);

      case

          s = (xi.*(.*eta + ).*(xi - ))/(*L2);

     case

         s = -((.*eta - ).*(xi.^ - ))/(*L2);

     case

         s = -(eta.*xi.*(xi + ))/L2;

     case

         s =-((*eta + ).*(xi.^ - ))/(*L2);

     case

         s = -(eta.*xi.*(xi - ))/L2;

     case

         s = (*eta.*(xi.^ - ))/L2;

      end

     end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 % subfunction phi

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

     function s = phi(xi,eta,i) % gradient of basis phi

     switch i

         case

             s = (eta.*xi.*(eta - ).*(xi - ))/;

         case

             s = (eta.*xi.*(eta - ).*(xi + ))/;

         case

             s = (eta.*xi.*(eta + ).*(xi + ))/;

         case

             s = (eta.*xi.*(eta + ).*(xi - ))/;

         case

             s = -(eta.*(xi.^ - ).*(eta - ))/;

         case

             s = -(xi.*(eta.^ - ).*(xi + ))/;

         case

             s = -(eta.*(xi.^ - ).*(eta + ))/;

         case

             s = -(xi.*(eta.^ - ).*(xi - ))/;

         case

             s = (eta.^ - ).*(xi.^ - );

     end

     end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 end

 function [elem2dof,edge,inDof] = c0dofP2(elem)

 totalEdge=sort([elem(:,[,]);elem(:,[,]);elem(:,[,]);elem(:,[,])],);

 [edge, i2, j] = myunique(totalEdge);

 N = max(elem(:));

 NT = size(elem,);

 NE = size(edge,);

 elem2edge = reshape(j,NT,);

 elem2dof = uint32([elem N+elem2edge (N+NE+:N+NE+NT)']);

 i1(j(*NT:-:)) = *NT:-:;

 i1 = i1';

 bdEdgeIdx = (i1 == i2);

 isBdDof = false(N+NE+NT,);

 isBdDof(edge(bdEdgeIdx,:)) = true;   % nodal

 idx = find(bdEdgeIdx);

 isBdDof(N+idx) = true;

 inDof = find(~isBdDof);

 end