how to find temperature of a wall in transient heat conduction in matlab

main.m
% Constants
k = 0.05; % Thermal conductivity [W/mK]
rho = 800; % Density [kg/m3]
Cp = 1800; % Heat capacity [J/kgK]

% Geometry
L = 1; % Length [m]
dx = 0.01; % Spacing [m]
N = L/dx + 1; % Number of nodes
x = linspace(0, L, N)'; % Node positions [m]

% Time
dt = 0.01; % Time step [s]
tFinal = 100; % Final time [s]
t = 0:dt:tFinal; % Time vector

% Initial condition
T0 = 20; % Initial temperature [C]
T = T0*ones(N, 1); % Temperature vector

% Boundary condition
T(N) = 100; % Constant temperature [C]

% Solver
r = k/(rho*Cp*dx^2); % Stability criterion
M = diag(2*r*ones(N-2, 1)) + diag(-r*ones(N-3, 1), 1) + diag(-r*ones(N-3, 1), -1); % Tridiagonal matrix
for i = 1:length(t)-1
    b = T(2:end-1) + r*(T(1:end-2) - 2*T(2:end-1) + T(3:end)); % RHS of the system
    T(2:end-1) = M\b; % Solve the system
end

% Plot
plot(x, T)
xlabel('Position [m]')
ylabel('Temperature [C]')
title('Temperature vs. Position')
913 chars
37 lines