% 2-D steady-state conduction
Tair=300; h=10; dx=0.25;k=1;Bi=h*dx/k;Ts=500;
figure(1); clf;

A=[-4   1  1  0  0  0  0  0;
    2  -4  0  1  0  0  0  0;
    1   0 -4  1  1  0  0  0;
    0   1  2 -4  0  1  0  0;
    0   0  1  0 -4  1  1  0;
    0   0  0  1  2 -4  0  1;
    0   0  0  0  2  0 -4-2*Bi 1;
    0   0  0  0  0  2  2 -4-2*Bi]
 
 bs=[-2; -1; -1; 0; -1; 0; -1; 0]
 bair=[0; 0; 0; 0; 0; 0; Tair; Tair]
 b=Ts*bs-2*Bi*bair
 
 T=inv(A)*b
 for i=1:5
    xx(i)=(i-1)*dx;
    for j=1:5
       yy(j)=(j-1)*dx;
       Z(i,j)=Ts;
    end
 end

 Z(1,2)=T(7);Z(1,3)=T(8);Z(1,4)=T(7);
 Z(2,2)=T(5);Z(2,3)=T(6);Z(2,4)=T(5);
 Z(3,2)=T(3);Z(3,3)=T(4);Z(3,4)=T(3);
 Z(4,2)=T(1);Z(4,3)=T(2);Z(4,4)=T(1);
 contourf(xx,yy,Z,10)
colormap('jet')
xlabel('x')
ylabel('y')


%  fig=figure(1);
%  print(fig,fullfile(pwd, 'Figure6_6.tif'),'-dtiff','-r600')