```clear;
format long;
fc=1;    % f constant
bc=0.;     % The Dirichlet boundary condition, constant
W=1;   % the width of the rectangle
L=2.0;    % the length of the rectangle
N=20;      % the number of elements in Y direction
M=40;       % the number of elements in X direction
dx = L/M;    % dx
dy = W/N;    % dy
D=(N+1)*(M+1);    % The dimension of the global stiffniss matrix or number of nodes
s=zeros(D,D);      % global stiffness matrix
v=zeros(N+1,M+1);    % the collocated results for plotting
xindex=0;
yindex=0;
x1=0; x2=0; x3=0; x4=0;   % (x1,x2,x3) and (x4,x2,x3) are summits of the two trangles in one rectangle element
j1=0; j2=0;
% generating the global stiffness matrix and load vectors
for yindex = 1:N
for xindex = 1:M
x1 = xindex + (yindex-1)*(M+1);
x2 = 1+ x1;
x3 = xindex + yindex*(M+1);
x4 = 1+ x3;

j1 = dx*dy;  % the Jacobian of the first trangle
s(x1,x1) = s(x1,x1)+(dx^2+dy^2)/(2*j1);
s(x1,x2) = s(x1,x2)-dy^2/(2*j1);
s(x1,x3) = s(x1,x3)-dx^2/(2*j1);
s(x2,x1) = s(x2,x1)-dy^2/(2*j1);
s(x2,x2) = s(x2,x2)+dy^2/(2*j1);
s(x2,x3) = s(x2,x3);
s(x3,x1) = s(x3,x1)-dx^2/(2*j1);
s(x3,x2) = s(x3,x2);
s(x3,x3) = s(x3,x3)+dx^2/(2*j1);
f(x1,1) = f(x1,1)+ fc*j1/6;
f(x2,1) = f(x2,1)+ fc*j1/6;
f(x3,1) = f(x3,1)+ fc*j1/6;   %  Process the first trangle in one block

j2 = dx*dy;   % the Jacobian of the second trangle
s(x4,x4) = s(x4,x4)+(dx^2+dy^2)/(2*j2);
s(x4,x2) = s(x4,x2)-dx^2/(2*j2);
s(x4,x3) = s(x4,x3)-dy^2/(2*j2);
s(x2,x4) = s(x2,x4)-dx^2/(2*j2);
s(x2,x2) = s(x2,x2)+dx^2/(2*j2);
s(x2,x3) = s(x2,x3);
s(x3,x4) = s(x3,x4)-dy^2/(2*j2);
s(x3,x2) = s(x3,x2);
s(x3,x3) = s(x3,x3)+dy^2/(2*j2);
f(x4,1) = f(x4,1)+ fc*j2/6;
f(x2,1) = f(x2,1)+ fc*j2/6;
f(x3,1) = f(x3,1)+ fc*j2/6;  %  Process the second trangle
end
end

% applying BC
for xindex=1:M+1
s(xindex,:) = zeros(1,D); s(xindex,xindex)=1; f(xindex,1)=bc;
s(xindex+N*(M+1),:) = zeros(1,D);
s(xindex+N*(M+1),xindex+N*(M+1))=1;
f(xindex+N*(M+1),1)=bc;
end         % the bottome and upper BC
for yindex=1:N-1
s(1+yindex*(M+1),:)=zeros(1,D);
s(1+yindex*(M+1),1+yindex*(M+1))=1;
f(1+yindex*(M+1),1)=bc;
s((1+yindex)*(M+1),:)=zeros(1,D);
s((1+yindex)*(M+1),(1+yindex)*(M+1))=1;
f((1+yindex)*(M+1),1)=bc;
end         % the left and right side BC

q=s\f;      % solve q

% generating visulized results
for yindex=1:N+1
for xindex=1:M+1
v(yindex,xindex) = q((M+1)*(yindex-1)+xindex);
end
end
[X,Y] = meshgrid(0:dx:L,0:dy:W);
[C,h] = contourf(X,Y,v);
clabel(C,h)
colormap cool
%axis equal tight
```