2d pde, implicit euler in time, central differencing in space. time loop with jacobi solver, discretized unknown ae aw an as ap pointers, to plot temperature on a plate a 0, 1, 2, 5, 10, 100 seconds in matlab

main.m
% Parameters
L = 1;          % Plate side length (m)
N = 20;         % Number of grid points in x and y directions
Tmax = 100;     % Maximum simulation time (s)
dt = 0.1;       % Time step size (s)
alpha = 1e-4;   % Thermal diffusivity (m^2/s)

% Calculated parameters
h = L/(N-1);    % Grid spacing (m)
x = linspace(0, L, N);   % Grid points x
y = linspace(0, L, N);   % Grid points y

% Initial conditions
T = zeros(N, N);   % Temperature matrix
T(:, 1) = 100;     % Left boundary condition
T(:, N) = 0;       % Right boundary condition
T(1, :) = 0;       % Top boundary condition
T(N, :) = 0;       % Bottom boundary condition

% Implicit Euler loop
for t = 0:dt:Tmax
    % Construct the coefficient matrix and RHS vector
    ae = -alpha*dt/(h^2)*ones(N-2, 1);
    aw = ae;
    an = ae;
    as = ae;
    ap = 1 + 4*alpha*dt/(h^2)*ones(N-2, N-2);
    b = T(2:N-1, 2:N-1);
    b(:, 1) = b(:, 1) + alpha*dt/(h^2)*T(2:N-1, 1);          % Neumann BC 1
    b(:, end) = b(:, end) + alpha*dt/(h^2)*T(2:N-1, end);      % Neumann BC 2
    b(1, :) = b(1, :) + alpha*dt/(h^2)*T(1, 2:N-1);            % Neumann BC 3
    b(end, :) = b(end, :) + alpha*dt/(h^2)*T(end, 2:N-1);      % Neumann BC 4
    
    % Jacobi Solver
    niter = 100;        % Number of Jacobi iterations
    for k = 1:niter
        Tnew = zeros(N, N);   % Initialize new temperature matrix
        Tnew(2:N-1, 2:N-1) = (-ae.*Tnew(2:N-1, 3:N) - aw.*Tnew(2:N-1, 1:N-2) ...
            - an.*Tnew(3:N, 2:N-1) - as.*Tnew(1:N-2, 2:N-1) + b)./ap;   % Update interior points
        T = Tnew;             % Update temperature matrix
    end
    
    % Plot
    contourf(x, y, T);
    colorbar;
    caxis([0 100]);
    xlabel('x (m)');
    ylabel('y (m)');
    title(sprintf('Temperature at t = %.1f s', t));
    drawnow;
end
1777 chars
52 lines