2d pde, implicit euler in time, central differencing in space. time loop with jacobi solver to plot temperature on a plate a 0, 1, 2, 5, 10, 100 seconds in matlab

main.m
% Plate size
Lx = 1;
Ly = 1;

% Number of spatial grid points
nx = 51;    % include boundary points
ny = 51;

dx = Lx/(nx-1);
dy = Ly/(ny-1);

% Spatial grid
x = linspace(0, Lx, nx);
y = linspace(0, Ly, ny);

[X, Y] = meshgrid(x, y);

% Thermal diffusivity
alpha = 0.01;

% Time variables
t0 = 0;
tend = 100;
nt = 10001;   % Number of time steps
dt = (tend - t0)/(nt-1);

% Initial conditions
T0 = 100 * sin(pi*X) .* sin(pi*Y);

% Boundary conditions
Tleft = 0;
Tright = 0;
Ttop = 0;
Tbottom = 0;

% Construct matrices for the implicit scheme
ex = ones(nx, 1);
ey = ones(ny, 1);
Ax = spdiags([ex, -2*ex, ex], [-1, 0, 1], nx, nx)/(dx^2);
Ay = spdiags([ey, -2*ey, ey], [-1, 0, 1], ny, ny)/(dy^2);
A = kron(Ay, speye(nx)) + kron(speye(ny), Ax);
B = speye(nx*ny) + alpha*dt*A;

% Time loop
T = T0(:);   % Initial temperature flattened to a column vector
for n = 1:nt
    % Solve the system of linear equations using Jacobi iteration
    for k = 1:20
        T = T + dt*alpha*A*T;
        T = jacobi(A, T, Tleft, Tright, Ttop, Tbottom, nx, ny);   % Jacobi iteration
    end
    % Plot temperature at specified times
    if any(n*dt == [0, 1, 2, 5, 10, 100])
        Tn = reshape(T, ny, nx);
        figure
        surf(X, Y, Tn)
        title(sprintf('Temperature at t = %.2f', n*dt))
        xlabel('x')
        ylabel('y')
        zlabel('T')
    end
end

function T = jacobi(A, T, Tleft, Tright, Ttop, Tbottom, nx, ny)
% Jacobi solver for Dirichlet boundary conditions
T = reshape(T, ny, nx);    % Reshape T to 2D grid
Told = T;

for i = 2:nx-1
    for j = 2:ny-1
        T(i, j) = (Told(i-1, j) + Told(i+1, j) + Told(i, j-1) + Told(i, j+1))/4;
    end
end

% Boundary conditions
T(1, :) = Ttop;
T(nx, :) = Tbottom;
T(:, 1) = Tleft;
T(:, ny) = Tright;

T = T(:);   % Flatten T to a column vector
end
1798 chars
83 lines