Please go through the attached code content on Finite difference method of steady state heat conduction problem and modify the MATLAB code (SteadyState. m) to simulate the temperature distribution at steady state

clear

clc

% close all

To=400; % Root temperature

TT=100; % ambient temperature

h=1000;

K=400;

DX=.1;

Bi=h*DX/K;% Biot Number

aa=-2*(Bi+1);

ab=-2*(Bi+2);

%% Matrix C: Righthand of the equation: [A]*[T]=[C]

cc=-TT*2*Bi;

C=[cc cc cc To cc 0 0 To cc 0 0 To cc cc cc To]';

%% Matrix A : Lefthand of the equation: [A]*[T]=[C]

A=zeros(16,16);

% Equation 1

A(1,1)=aa;

A(1,5)=1;

A(1,2)=1;

% Equation 2

A(2,2)=ab;

A(2,6)=2;

A(2,1)=1;

A(2,3)=1;

% Equation 3

A(3,3)=ab;

A(3,7)=2;

A(3,2)=1;

A(3,4)=1;

% Equation 4

A(4,4)=1;

% Equation 5

A(5,5)=ab;

A(5,6)=2;

A(5,9)=1;

A(5,1)=1;

% Equation 6

A(6,6)=-4;

A(6,5)=1;

A(6,10)=1;

A(6,2)=1;

A(6,7)=1;

% Equation 7

A(7,7)=-4;

A(7,6)=1;

A(7,11)=1;

A(7,3)=1;

A(7,8)=1;

% Equation 8

A(8,8)=1;

% Equation 9

A(9,9)=ab;

A(9,10)=2;

A(9,13)=1;

A(9,5)=1;

% Equation 10

A(10,10)=-4;

A(10,14)=1;

A(10,9)=1;

A(10,6)=1;

A(10,11)=1;

% Equation 11

A(11,11)=-4;

A(11,15)=1;

A(11,10)=1;

A(11,7)=1;

A(11,12)=1;

% Equation 12

A(12,12)=1;

% Equation 13

A(13,13)=aa;

A(13,14)=1;

A(13,9)=1;

% Equation 14

A(14,14)=ab;

A(14,10)=2;

A(14,13)=1;

A(14,15)=1;

% Equation 15

A(15,15)=ab;

A(15,11)=2;

A(15,14)=1;

A(15,16)=1;

% Equation 16

A(16,16)=1;

%% Solve for Temperature

T=linsolve(A,C);

TR=reshape(T,[4,4]);

TR=TR';

X=(0:3)*DX;

Y=X;

[XX,YY]=meshgrid(X,Y);

% plot temperature distribution

figure(2)

contourf(XX,YY,TR)

colormap(jet)

xlabel('X')

xlabel('Y')

colorbar

caxis([TT To])