a=0;b=1;
c=0;d=1;
n=32;m=32;
h=(b-a)/n;
l=(d-c)/m;
NK=(n+1)*(m+1);%节点的总个数%
NE=2*n*m;%剖分的小三角形的总个数%
x=zeros(n+1,1);
y=zeros(n+1,1);
for r=1:n+1
    x(r)=a+h*(r-1);
end
for g=1:m+1
    y(g)=c+l*(g-1);
end%形成坐标轴上的坐标%
A=zeros(NK,2);%节点、节点横坐标、节点纵坐标所对应的总矩阵%
X=zeros(NK,1);
Y=zeros(NK,1);
for k=1:NK
    A(k)=k;
    q=fix(k/(n+1));%最靠近0的整数%
    s=rem(k,n+1);%余数%
    if (s==0&&k<=NK-1)%可以整除，此时节点在最后一列%
        X(k)=b;
        Y(k)=y(q);
    elseif (s~=0&&k<=NK-1);
        X(k)=x(s);%不能整除，此时节点的位置判断如右%
        Y(k)=y(q+1);
    else X(k)=b;Y(k)=d;    
    end
    A(k,1)=X(k);
    A(k,2)=Y(k);
end
B=zeros(NE,3);%三角形编号、三角形三个点分别对应的节点编号的总矩阵%
for p=1:NE%三角形的三个顶点所对应的整体的节点编码表示出来%
   Q1=fix(p/(2*n));%hangshu最靠近0的整数，三角形所在的行%
   S1=rem(p,2*n);%lieshu余数%
   S2=rem(p,2);
   if (S1==0)%每一行的最后一个三角形
        B(p,2)=Q1*(n+1); B(p,1)= B(p,2)-1;B(p,3)= B(p,2)+n+1;
   else
       if (S2==0)
            B(p,1)=Q1*(n+1)+S1/2; B(p,2)= B(p,1)+1;B(p,3)= B(p,2)+n+1;
       else
            B(p,1)=Q1*(n+1)+(S1+1)/2;B(p,3)=B(p,1)+n+1;B(p,2)=B(p,3)+1;
       end
   end
 end
   %构造逼近函数%
MP=zeros(3,3);
%DZ2=zeros(NK,NK);
DZ=zeros(NK,NK);
S2=zeros(NK,1);
for p=1:NE
   S1=zeros(NK,1);
   DZ1=zeros(NK,NK);
   m1= B(p,1);m2= B(p,2);m3= B(p,3);
   MP=[A(m1,1) A(m1,2) 1;A(m2,1) A(m2,2) 1;A(m3,1) A(m3,2) 1]; 
   sp=det(MP);
   T=magic(3);
   H=[MP(1,1) MP(1,2);
    MP(2,1) MP(2,2);
    MP(3,1) MP(3,2);
    MP(1,1) MP(1,2);
    MP(2,1) MP(2,2);];
for i=1:3
    T(i,1)=H(i+1,2)-H(i+2,2);
    T(i,2)=-H(i+1,1)+H(i+2,1);
    T(i,3)=H(i+1,1)*H(i+2,2)-H(i+2,1)*H(i+1,2);
end
   D=magic(3);  
     N=1/sp*T;%第p个三角形三个顶点的基函数系数矩阵%
        for i=1:3
            for j=1:3
                D(i,j)=0.5*sp*(N(i,1)*N(j,1)+N(i,2)*N(j,2)) ; %左边积分值 
            end
        end
     for  i=1:3
         for j=1:3
     DZ1(B(p,i),B(p,j))=D(i,j);
         end
        S1(B(p,i))=(1/3)*sp;
     end
     DZ= DZ+DZ1;%总刚度矩阵
     S2=S1+S2;
     %右端的积分值放在矩阵中，用于下面划去边界值所对应的行
end
DZ=[DZ S2];
dl=NK-2*(n+m);
j1=1;j2=1;j3=1;j4=1;j=1;
u1=zeros(n,1);
u2=zeros(m+1,1);
u3=zeros(m,1);
u4=zeros(n-1,1);
DZZ=zeros(dl,NK+1);
for i=1:NK     %划掉行
     aa=rem(i,n+1);
     bb=fix(i/(n+1));
     if  (i<=n)
     u1(j1)=A(i,1)^2;
        j1=j1+1;
     elseif (aa==0)
       u2(j2)=b.^2;
       j2=j2+1;
      elseif (i~=1&&aa==1)
           u3(j3)=a^2;
           j3=j3+1;
      elseif (i~=NK&&i~=NK-n&&bb==m)
           u4(j4)=A(i,1)^2;
           j4=j4+1;
      else 
          DZZ(j,:)=DZ(i,:);%magic
          j=j+1;
     end
end
 S=DZZ(:,NK+1);%形成右端数
 DZL=zeros(dl,dl);
 F1=zeros(dl,n);
 F2=zeros(dl,m+1);
 F3=zeros(dl,m);
 F4=zeros(dl,n-1);
 j1=1;j2=1;j3=1;j4=1;j=1;
 for i=1:NK    %划掉列
     aa=rem(i,n+1);
     bb=fix(i/(n+1));
     if  (i<=n)
      F1(:,j1)=DZZ(:,i);
      j1=j1+1;
     elseif (aa==0)
      F2(:,j2)=DZZ(:,i);
          j2=j2+1;
      elseif (i~=1&&aa==1)
          F3(:,j3)=DZZ(:,i);
           j3=j3+1;
      elseif (i~=NK&&i~=NK-n&&bb==m)
           F4(:,j4)=DZZ(:,i);
           j4=j4+1;
      else 
          DZL(:,j)=DZZ(:,i);%magic
          j=j+1;
     end
  end
U=[u1;u2;u3;u4];
F=[F1,F2,F3,F4];
SL=-S-F*U;
UL=pinv(DZL)*SL;
ur=zeros(dl,1);
for i=1:dl
    ar=rem(i,n-1);
    if(ar~=0)
        ur(i)=x(ar+1).^2;
    else
        ur(i)=x(n).^2;
    end
end
err=(abs(UL-ur))./ur;
error=max(err);
UH=UL'*DZL*UL;
%DZL*ur+F*U