I have a problem with 2D heat conduction develop a computer code to solve the 2D heat conduction equation:∂2T /∂x2 + ∂2T/ ∂y2 = 0, gradient don't change when iterations increase, please help me

clc; clear; close all;

% === COMMON DOMAIN SETUP ===
L = 1.0;
Nx = 41; Ny = 41;
dx = L / (Nx - 1);
dy = dx;
x = linspace(0, L, Nx);
y = linspace(0, L, Ny);
[X, Y] = meshgrid(x, y);
alpha = 1.0;
tol = 1e-5;

% === TRACK ITERATIONS ===
iterations_to_plot = [1 2 3 5 10 20 40 80 120 160 200];

%% === 1. STEADY-STATE SOLVER (Gauss-Seidel) ===

T_steady = zeros(Ny, Nx);
    % Apply boundary conditions
    T_steady(:,1)   = 0;  % Left
    T_steady(:,end) = 0;  % Right
    T_steady(1,:)   = 0;  % Bottom
    T_steady(end,:) = 1;  % Top

residual = Inf;
iter = 0;
max_iter = 10000;
grad_vec_steady = [];  % Khởi tạo mảng rỗng để lưu max gradient
T_steady_snapshots = cell(length(iterations_to_plot), 1);
snap_idx = 1;

while residual > tol && iter < max_iter
    T_old = T_steady;
    
    % Gauss-Seidel update
    for j = 2:Ny-1
        for i = 2:Nx-1
            T_steady(j,i) = 0.25 * ( ...
                T_steady(j+1,i) + T_steady(j-1,i) + ...
                T_steady(j,i+1) + T_steady(j,i-1));
        end
    end

    residual = max(max(abs(T_steady - T_old)));
    iter = iter + 1;

    % % --- Compute gradient magnitude ---
    [Tx, Ty] = gradient(T_steady, dx, dy);
    grad_mag = sqrt(Tx.^2 + Ty.^2);
    grad_vec_steady(iter) = max(grad_mag(:));

    % --- Save snapshot if needed ---
    if snap_idx <= length(iterations_to_plot) && iter == iterations_to_plot(snap_idx)
        T_steady_snapshots{snap_idx} = T_steady;
        snap_idx = snap_idx + 1;
    end
end

fprintf(' Steady-state converged in %d iterations. Final residual: %.2e\n', iter, residual);

 Steady-state converged in 905 iterations. Final residual: 9.97e-06

% --- Plot Max Temperature Gradient vs Iteration ---
figure;
plot(1:iter, grad_vec_steady(1:iter), 'LineWidth', 1.5);
xlabel('Iteration');
ylabel('Max Temperature Gradient');
title('Maximum Temperature Gradient vs Iteration (Steady-State)');
grid on;

% --- Plot: Converged Steady-State Temperature Distribution ---
figure;
contourf(X, Y, T_steady, 20, 'LineColor', 'none');
colorbar;
xlabel('x');
ylabel('y');
title('Converged Steady-State Temperature Distribution (Gauss-Seidel)');

Answer 1

[Torsten]

The maximum gradient is always in the left and right upper corner point with value (1-0)/(1/40) = 40. So it won't change in the course of the iteration.