function [A,b]=fem(u,v,a0,b0)
% a 是x軸上的邊界長度 ,b是y軸上的邊界長度
% m是x軸上的等分數，n是y軸上的等分數，Nnode為節點總數
% node全局編碼數組，為一維數組，存放節點全局編號
% Elematrix單元號與節點局部編碼和全局編碼的二維數組
% clear all
% clc

% u=32;v=32;a0=1;b0=1;
hx=a0/u;hy=b0/v;
N_node=(u+1)*(v+1); N_elem=2*u*v;
delta_e=hx*hy*0.5;

%生成節點總體編碼與各單元的節點局部編碼的矩陣
%本程序中的節點編碼及其排序是根據《偏微分方程數值解法》（陸金甫，關治）P240中例子 

for e=1:N_elem
    alphax(e)=1;
    alphay(e)=1;
    beta(e)=0;
    t=2*u;
    if mod(e,t)==0
        t=e/t-1;
        x(e)=ceil(e/2)+t;
        M_ju_zong(1,e)=x(e)+1;
        M_ju_zong(2,e)=x(e)+u+2;
        M_ju_zong(3,e)=x(e);
    else
        x(e)=ceil(e/2)+floor(e/t);
        if mod(e,2)==1
            M_ju_zong(1,e)=x(e)+u+1;
            M_ju_zong(2,e)=x(e);
            M_ju_zong(3,e)=x(e)+u+2;
        else
            M_ju_zong(1,e)=x(e)+1;
            M_ju_zong(2,e)=x(e)+u+2;
            M_ju_zong(3,e)=x(e);
        end
    end
end


%初始化總體剛度矩陣
for i=1:N_node
    b(i)=0;
    for j=1:N_node
        A(i,j)=0;
    end
end

%計算ai,aj,am,bi,bj,bm

for e=1:N_elem
      i(e)=M_ju_zong(1,e);
      j(e)=M_ju_zong(2,e);
      m(e)=M_ju_zong(3,e);
%確定i（e）的坐標
      r=ceil(i(e)/(u+1));%第i個節點的行號(取右邊整數)
      s=mod(i(e),u+1);   %第i個節點的列號
      if s==0
         s=u+1;
      end
      node(i(e))=(u+1)*(r-1)+s;
      x(i(e))=(s-1)*hx;
      y(i(e))=(r-1)*hy;
%       fe(i(e))=0;
      fe(i(e))=-2*pi^2*sin(pi*x(i(e)))*sin(pi*y(i(e)));
%確定j（e）的坐標
      r=ceil(j(e)/(u+1));%第i個節點的行號(取右邊整數)
      s=mod(j(e),u+1);   %第i個節點的列號
      if s==0
         s=u+1;
      end
      node(j(e))=(u+1)*(r-1)+s;
      x(j(e))=(s-1)*hx;
      y(j(e))=(r-1)*hy;
%       fe(i(e))=0;
      fe(j(e))=-2*pi^2*sin(pi*x(j(e)))*sin(pi*y(j(e)));
%確定m（e）的坐標
      r=ceil(m(e)/(u+1));%第i個節點的行號(取右邊整數)
      s=mod(m(e),u+1);   %第i個節點的列號
      if s==0
         s=u+1;
      end
      node(m(e))=(u+1)*(r-1)+s;
      x(m(e))=(s-1)*hx;
      y(m(e))=(r-1)*hy;
%       fe(i(e))=0;
      fe(m(e))=-2*pi^2*sin(pi*x(m(e)))*sin(pi*y(m(e)));
      
     be(1)=y(j(e))-y(m(e));
     be(2)=y(m(e))-y(i(e));
     be(3)=y(i(e))-y(j(e));
     ce(1)=x(m(e))-x(j(e));
     ce(2)=x(i(e))-x(m(e));
     ce(3)=x(j(e))-x(i(e));
end
deltae=0.5*(be(1)*ce(2)-be(2)*ce(1));
%計算delta_ij

for i=1:3
    for j=1:3
        if (i~=j)
            del_ij=1.0;
        else
            del_ij=0.0;
        end
    end
end
    
for e=1:N_elem
%計算單元剛度矩陣
     for i=1:3
         b(i)=deltae*fe(m(e))/3;
         for j=1:3
             kk(i,j)=(alphax(e)*be(i)*be(j)+alphay(e)*ce(i)*ce(j))/(4*deltae)+beta(e)*(1.0+del_ij)*deltae/12;
         end
     end
%合成單元剛度矩陣到剛度矩陣
     for i=1:3
        for j=1:3
            A(M_ju_zong(i,e),M_ju_zong(j,e))=A(M_ju_zong(i,e),M_ju_zong(j,e))+kk(i,j);
            b(M_ju_zong(i,e))=b(M_ju_zong(i,e))+b(i);
        end
     end
end

%邊界條件的處理

 A(1:u+1,:)=[];%去掉矩陣前u+1行
 A(:,1:u+1)=[];%去掉矩陣前u+1列
 A=A(1:N_node-2*(u+1),:);%去掉矩陣后u+1行
 A=A(:,1:N_node-2*(u+1));%去掉矩陣后u+1列
 b(1:u+1)=[];%去掉右端向量前u+1行
 b(N_node-2*(u+1)+1:N_node-(u+1))=[];%去掉右端向量后u+1行

 %去掉左邊界結點所對應的行列
 
 g=N_node-2*(u+1);p=0;
 for i=1:(g-p)
     if mod(i+p,u+1)==1&i+p<=g
         A(i,:)=[];
         A(:,i)=[];
         b(i)=[];
         p=p+1;
     end
 end
%去掉右邊界結點所對應的行列

g1=g-v+1;p=0;
 for i=1:(g1-p)
     if mod(i+p,u)==0&i+p<=g1
         A(i,:)=[];
         A(:,i)=[];
         b(i)=[];
         p=p+1;
     end
 end
 b=b';
% size(k);size(b);
% spy(k);
% disp(b');
% k;