Here is a sample code in MATLAB that solves the 2D heat equation with a Jacobi solver and implicit Euler time scheme, using central differencing in space:
% Parameters
diffusivity = 0.01;
final_time = 100;
delta_t = 0.5;
delta_x = 0.05;
delta_y = 0.05;
x = 0:delta_x:1;
y = 0:delta_y:1;
Nx = numel(x);
Ny = numel(y);
Nt = floor(final_time / delta_t) + 1;
T = zeros(Nx, Ny, Nt);
% Initial condition
T(:, :, 1) = ones(Nx, Ny);
% Boundary conditions
T(1, :, :) = 0;
T(Nx, :, :) = 0;
T(:, 1, :) = 0;
T(:, Ny, :) = 0;
% Jacobi solver
for n = 2:Nt
T(:, :, n) = T(:, :, n-1);
for i = 2:Nx-1
for j = 2:Ny-1
T(i, j, n) = ( T(i+1, j, n-1) + T(i-1, j, n-1) ...
+ T(i, j+1, n-1) + T(i, j-1, n-1) ) / 4 ...
- diffusivity * delta_t / delta_x^2 ...
* ( T(i+1, j, n-1) - 2*T(i, j, n-1) + T(i-1, j, n-1) ...
+ T(i, j+1, n-1) - 2*T(i, j, n-1) + T(i, j-1, n-1) );
end
end
end
% Plotting
times = [0 1 2 5 10 100];
figure;
for i = 1:numel(times)
[~, t] = min(abs(delta_t*(1:Nt) - times(i)));
subplot(2, 3, i);
surf(x, y, T(:, :, t)');
shading interp;
title(['t = ' num2str(delta_t * t)]);
end
Note that the Jacobi solver only updates the temperature field after a complete iteration over all grid points. This can be quite slow, especially for large systems or long simulation times. More advanced solvers like Gauss-Seidel or successive over-relaxation (SOR) can be used to speed up the computations at the cost of increased complexity.
The code above generates a 2D plot of the temperature field at selected time points. The surf function allows for various visualization options, such as changing the colormap, adding a colorbar, or rotating the view.