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

main.m
% Define plate size and time step
Lx = 1; % plate length in x-direction
Ly = 1; % plate length in y-direction
Nx = 10; % number of grid points in x-direction
Ny = 10; % number of grid points in y-direction
dx = Lx/(Nx-1); % grid spacing in x-direction
dy = Ly/(Ny-1); % grid spacing in y-direction
dt = 0.01; % time step

% Define initial temperature distribution
T0 = zeros(Ny,Nx); % initialize temperature matrix
T0(1,:) = 100; % bottom boundary
T0(Ny,:) = 0; % top boundary
T0(:,1) = 0; % left boundary
T0(:,Nx) = 0; % right boundary

% Define diffusion coefficient
D = 1; % thermal diffusivity

% Define coefficient matrices
Ax = toeplitz([2 -1 zeros(1,Nx-2)])/(dx^2); % x-direction matrix
Ay = toeplitz([2 -1 zeros(1,Ny-2)])/(dy^2); % y-direction matrix
A = kron(Ay,speye(Nx)) + kron(speye(Ny),Ax); % 2D matrix
I = speye(Nx*Ny); % identity matrix

% Define implicit scheme
Tn = T0(:); % reshape initial temperature matrix into a column vector
for t = 0:dt:100
    b = Tn; % RHS of system of equations
    Tn = jacobi_solve(I - (dt*D/dx^2)*A,b); % solve system of equations
    if mod(t,1)==0 % plot temperature every specified time step
        T = reshape(Tn,Ny,Nx); % reshape temperature column vector into a matrix
        surf(linspace(0,Lx,Nx),linspace(0,Ly,Ny),T) % plot temperature surface
        title(['Temperature at t=',num2str(t),' seconds']) % title contains time
        xlabel('x') % x-axis label
        ylabel('y') % y-axis label
        zlabel('Temperature') % z-axis label
        drawnow % update plot
    end
end
1540 chars
41 lines

main.m
function x = jacobi_solve(A,b,max_iter,tol)
% Solves the system of equations Ax=b using Jacobi method
% Inputs:
%   A: coefficient matrix
%   b: right-hand side vector
%   max_iter (optional): maximum number of iterations (default=1000)
%   tol (optional): convergence tolerance (default=1e-6)
% Output:
%   x: solution vector

% Set default values for optional inputs
if nargin < 3 || isempty(max_iter)
    max_iter = 1000;
end
if nargin < 4 || isempty(tol)
    tol = 1e-6;
end

% Initialize variables
n = length(b);
x = zeros(n,1);
x_old = x;

% Jacobi iteration
for k = 1:max_iter
    for i = 1:n
        x(i) = (b(i) - A(i,:) * x_old) / A(i,i);
    end
    if norm(x - x_old) / norm(x_old) < tol % check for convergence
        break
    end
    x_old = x;
end

% Display warning if maximum iterations is reached and no convergence
if k == max_iter
    warning('Maximum iterations reached without convergence')
end
end
923 chars
40 lines