Errors in transient heat transfer node

% Simulation Parameters
dt =10; % Time step [s]
t_max = 2000; % Maximum simulation time [s]
time_steps = t_max / dt;
Nx = 34; Ny = 40; % Grid dimensions
dx = 1; dy = 1; % Uniform grid spacing
h = ones(Ny, Nx) * 15;
h_air=15;
h_water=300;
h_in=2;
% Initialize material map and properties
material_map = 10 * ones(Ny, Nx); % Initialize material map
  % ABS
material_map(1:13, 6:8) = 1;
material_map(1:13, 28:30) = 1;
material_map(1:4, 9:27) = 1;
% Handle
material_map(29:30, 6:9) = 2;
material_map(20:28, 6:7) = 2;
% Vacuum
material_map(10:31,11) = 3; % Row 11, Columns 10 to 31
material_map(10, 11:12) = 3; % Row 10, Columns 11 to 12
material_map(6:10, 12) = 3;  % Rows 6 to 10, Column 12
material_map(6, 12:24) = 3;  % Row 6, Columns 12 to 24
material_map(6:10, 24) = 3;  % Rows 6 to 10, Column 24
material_map(10, 24:25) = 3; % Row 10, Columns 24 to 25
material_map(10:31,25) = 3;
% Water
material_map(12:30,13) = 4; % Row 13, Columns 12 to 30
material_map(8:30, 14) = 4; % Rows 8 to 30, Column 14
material_map(8:26, 15:21) = 4; % Rows 8 to 26, Columns 15 to 21
material_map(8:30, 22) = 4; % Rows 8 to 30, Column 22
material_map(12:30, 23) = 4; % Rows 12 to 30, Column 23
% Ice
material_map(27:30,15:21 ) = 5; % Rows 15 to 21, Columns 27 to 30
% Internal Air
material_map(31:34, 15:21) = 6; % Rows 31 to 34, Columns 15 to 21
material_map(31:32, 14) = 6;    % Rows 31 to 32, Column 14
material_map(31, 13) = 6;       % Row 31, Column 13
material_map(31, 23) = 6;       % Row 31, Column 23
material_map(31:32, 22) = 6;    % Rows 31 to 32, Column 22
% Lid
material_map(35:36, 14:22) = 7; % Rows 35 to 36, Columns 14 to 22
% Stainless Steel (Inside)
material_map(33, 13:14) = 8;    % Row 33, Columns 13 to 14
material_map(33, 22:23) = 8;    % Row 33, Columns 22 to 23
material_map(32, 23:24) = 8;    % Row 32, Columns 23 to 24
material_map(32, 12:13) = 8;    % Row 32, Columns 12 to 13
material_map(11:31, 12) = 8;    % Rows 11 to 31, Column 12
material_map(11:31, 24) = 8;    % Rows 11 to 31, Column 24
material_map(11, 23:24) = 8;    % Row 11, Columns 23 to 24
material_map(11, 12:13) = 8;    % Row 11, Columns 12 to 13
material_map(7:10, 13) = 8;     % Rows 7 to 10, Column 13
material_map(7:10, 23) = 8;     % Rows 7 to 10, Column 23
material_map(7, 13:23) = 8;     % Row 7, Columns 13 to 23
% Stainless Steel (Outside) 위치 정의
material_map(34, 12:14) = 9;    % Row 34, Columns 12 to 14
material_map(34, 22:24) = 9;    % Row 34, Columns 22 to 24
material_map(33, 11:12) = 9;    % Row 33, Columns 11 to 12
material_map(33, 24:25) = 9;    % Row 33, Columns 24 to 25
material_map(32, 25:26) = 9;    % Row 32, Columns 25 to 26
material_map(32, 10:11) = 9;    % Row 32, Columns 10 to 11
material_map(9:32, 26) = 9;     % Rows 9 to 32, Column 26
material_map(31:32, 11) = 9;    % Rows 31 to 32, Column 11
material_map(9:32, 10) = 9;     % Rows 9 to 32, Column 10
material_map(5:9, 11) = 9;      % Rows 5 to 9, Column 11
material_map(5:9, 25) = 9;      % Rows 5 to 9, Column 25
material_map(5, 11:25) = 9;     % Row 5, Columns 11 to 25
k(material_map == 1) = 0.23; % ABS
k(material_map == 2) = 0.3; % Handle
k(material_map == 3) = 0.0003; % Vacuum
k(material_map == 4) = 0.575; % Water
k(material_map == 5) = 2.2; % Ice
k(material_map == 6) = 0.025; % Internal Air
k(material_map == 7) = 0.2; % Lid
k(material_map == 8) = 16.2; % Stainless Steel (Inside)
k(material_map == 9) = 15.2; % Stainless Steel (Outside)
k(material_map == 10) = 0.025; % External Air
rho(material_map == 1) = 1040;
rho(material_map == 2) = 950;
rho(material_map == 3) = 0; % Vacuum
rho(material_map == 4) = 1000;
rho(material_map == 5) = 917;
rho(material_map == 6) = 1.2;
rho(material_map == 7) = 950;
rho(material_map == 8) = 8000;
rho(material_map == 9) = 8000;
rho(material_map == 10) = 1.2;
c(material_map == 1) = 1.5;
c(material_map == 2) = 2.0;
c(material_map == 3) = 0; % Vacuum
c(material_map == 4) = 4.186;
c(material_map == 5) = 2.06;
c(material_map == 6) = 1.005;
c(material_map == 7) = 1.75;
c(material_map == 8) = 0.5;
c(material_map == 9) = 0.5;
c(material_map == 10) = 1.005;
T(material_map == 1) = 21;
T(material_map == 2) = 19;
T(material_map == 3) = 14; % Vacuum
T(material_map == 4) = 4;
T(material_map == 5) = -15;
T(material_map == 6) = 8;
T(material_map == 7) = 15;
T(material_map == 8) = 5;
T(material_map == 9) = 19;
T(material_map == 10) = 20; % External
% Initialize material properties
[k, rho, c, T] = initialize_material_properties(material_map, Nx, Ny);
% Constants
T_infinite = 20; % Ambient temperature
sigma = 5.67e-8; % Stefan-Boltzmann constant
epsilon = 0.9; % Emissivity
video = VideoWriter('Temperature_Simulation.avi'); % 동영상 파일 이름
video.FrameRate = 10; % 프레임 속도 설정
open(video); % 동영상 파일 열기
% Simulation loop
all_temperatures = zeros(Ny, Nx, t_max / dt); % Ny x Nx 크기에 TimeSteps 추가
for t = 1:t_max
    % 새로운 온도 저장 배열
    T_new = T;
    % 내부 노드 및 경계 처리
    for i = 2:Ny-1
        for j = 2:Nx-1
             if rho(i, j) > 0 && c(i, j) > 0
            % 물질 내부의 전도 (Fourier's Law)
            Txx = (T(i+1, j) - 2*T(i, j) + T(i-1, j)) / dx^2;
            Tyy = (T(i, j+1) - 2*T(i, j) + T(i, j-1)) / dy^2;
            alpha = k(i, j) / (rho(i, j) * c(i, j)); % 열확산계수
            T_cond_internal = alpha * (Txx + Tyy);
            % 물질 간 경계에서의 열 전달
            T_boundary = 0;
            T_new(i, j) = T(i, j) + dt * T_cond_internal;
            end
           % 1. ABS - External Air (대류)
            if material_map(i, j) == 1 && is_material_boundary(i, j, material_map, 10)
                T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
                T_new(i, j) = T(i, j) + dt * T_conv;
            end
            % 2. ABS - External Stainless Steel (전도)
            if material_map(i, j) == 1 && is_material_boundary(i, j, material_map, 9)
                k_eff = 2 * k(i, j) * k(i+1, j) / (k(i, j) + k(i+1, j)); % Effective conductivity
                T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
                T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
            end
            % 3. External Stainless Steel - External Air (대류)
            if material_map(i, j) == 9 && is_material_boundary(i, j, material_map, 10)
                T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
                T_new(i, j) = T(i, j) + dt * T_conv;
            end
            % 4. External Stainless Steel - Vacuum (전도)
            if material_map(i, j) == 9 && is_material_boundary(i, j, material_map, 3)
                k_eff = k(i+1, j); % Vacuum's conductivity
                T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
                T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
            end
            % 5. Internal Stainless Steel - Water (대류)
            if material_map(i, j) == 8 && is_material_boundary(i, j, material_map, 4)
                T_conv = h_water * (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
                T_new(i, j) = T(i, j) + dt * T_conv;
            end
            % 6. Internal Stainless Steel - Vacuum (전도)
            if material_map(i, j) == 8 && is_material_boundary(i, j, material_map, 3)
                k_eff = k(i+1, j); % Vacuum's conductivity
                T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
                T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
            end
            % 7. Water - Internal Air (대류)
            if material_map(i, j) == 4 && is_material_boundary(i, j, material_map, 6)
                T_conv = 50 * (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
                T_new(i, j) = T(i, j) + dt * T_conv;
            end
            % 8. Water - Ice (전도)
            if material_map(i, j) == 4 && is_material_boundary(i, j, material_map, 5)
                k_eff = 2 * k(i, j) * k(i+1, j) / (k(i, j) + k(i+1, j));
                T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
                T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
            end
            % 9. Ice - Internal Air (대류)
            if material_map(i, j) == 5 && is_material_boundary(i, j, material_map, 6)
                T_conv = h_in * (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
                T_new(i, j) = T(i, j) + dt * T_conv;
            end
            % 10. Internal Air - Lid (대류)
            if material_map(i, j) == 6 && is_material_boundary(i, j, material_map, 7)
                T_conv = h_in* (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
                T_new(i, j) = T(i, j) + dt * T_conv;
            end
            % 11. Lid - External Stainless Steel (전도)
            if material_map(i, j) == 7 && is_material_boundary(i, j, material_map, 9)
                k_eff = 2 * k(i, j) * k(i+1, j) / (k(i, j) + k(i+1, j));
                T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
                T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
            end
            % 12. Lid - External Air (대류)
            if material_map(i, j) == 7 && is_material_boundary(i, j, material_map, 10)
                T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
                T_new(i, j) = T(i, j) + dt * T_conv;
            end
            % 13. Handle - External Stainless Steel (전도)
            if material_map(i, j) == 2 && is_material_boundary(i, j, material_map, 9)
                k_eff = 2 * k(i, j) * k(i+1, j) / (k(i, j) + k(i+1, j));
                T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
                T_new(i, j) = T(i, j) + dt * T_cond / (rho(i, j) * c(i, j));
            end
            % 14. Handle - External Air (대류)
            if material_map(i, j) == 2 && is_material_boundary(i, j, material_map, 10)
                T_conv = h(i, j) * (T_infinite - T(i, j)) / (rho(i, j) * c(i, j));
                T_new(i, j) = T(i, j) + dt * T_conv;
            end 
             if material_map(i, j) == 4 % 물
            % 상하좌우 인접 노드와의 대류 효과
            T_up = T(i-1, j);
            T_down = T(i+1, j);
            T_left = T(i, j-1);
            T_right = T(i, j+1);
            if is_material_boundary(i, j, material_map, 6)
    if material_map(i, j) == 4  % 물과 공기
        T_conv = h_water * (T(i+1, j) - T(i, j)) / (rho(i, j) * c(i, j));
    else % 기본 전도
        T_cond = k_eff * (T(i+1, j) - T(i, j)) / dx^2;
    end
    T_new(i, j) = T(i, j) + dt * (T_cond + T_conv);
end
            % 대류에 의한 열 이동 추가
            T_conv_internal = h_water/ (rho(i, j) * c(i, j)) * ...
                              ((T_up - T(i, j)) + (T_down - T(i, j)) + ...
                               (T_left - T(i, j)) + (T_right - T(i, j)));
            T_new(i, j) = T(i, j) + dt * T_conv_internal;
        end
            % 내부 전도와 경계에서의 열 전달을 모두 반영하여 온도 계산
            T_new(i, j) = T(i, j) + dt * (T_cond_internal + T_boundary);

        end
    end
    % 추가 경계 조건 적용
    % 하단 경계 (대류 경계, T_infinite = 20°C)
    for j = 1:Nx
        T_conv_bottom = h_air * (T_infinite - T(Ny, j)) / (rho(Ny, j) * c(Ny, j));
        T_new(Ny, j) = T(Ny, j) + dt * T_conv_bottom;
    end
    % 왼쪽 경계 (대류 경계, T_infinite = 20°C)
    for i = 1:Ny
        T_conv_left = h_air * (T_infinite - T(i, 1)) / (rho(i, 1) * c(i, 1));
        T_new(i, 1) = T(i, 1) + dt * T_conv_left;
    end
    % 상단 경계 (단열 경계)
    T_new(1, :) = T_new(2, :);
    % 우측 경계 (고정 온도)
    T_new(:, Nx) = 400;
    % 온도 업데이트
    T = T_new;
  if mod(t, 10) == 0 % High-resolution colormap
    imagesc(flipud(T), [-20, 40]); % Show the initial temperature field
    colorbar;
    title(['Time = ', num2str(t*dt), ' seconds']);
    xlabel('X-Axis (Nodes)');
    ylabel('Y-Axis (Nodes)');
    drawnow; % Ensure the figure updates properly
    frame = getframe(gcf); % Capture the current figure
    writeVideo(video, frame); % Save the frame to the video
    if mod(t, 500) == 0
            writematrix(T, ['Node_Temperatures_Time_', num2str(t * dt), '.csv']);
            disp(['Saved node temperatures at Time = ', num2str(t * dt), ' seconds']);
        end
  end
  if t * dt >= t_max
        break; % 종료
    end
end
writematrix(T, 'Final_Node_Temperatures.csv');
disp('Final node temperatures saved to Final_Node_Temperatures.csv');

    close(video);
% Local Function: is_material_boundary
function is_boundary = is_material_boundary(i, j, material_map, target_material)
    % 기본적으로 경계 아님
    is_boundary = false;
    % 인접 노드를 검사 (상하좌우)
    if i > 1 && material_map(i-1, j) == target_material % 위쪽
        is_boundary = true;
    elseif i < size(material_map, 1) && material_map(i+1, j) == target_material % 아래쪽
        is_boundary = true;
    elseif j > 1 && material_map(i, j-1) == target_material % 왼쪽
        is_boundary = true;
    elseif j < size(material_map, 2) && material_map(i, j+1) == target_material % 오른쪽
        is_boundary = true;
    end
end
% Local Function: initialize_material_properties
function [k, rho, c, T] = initialize_material_properties(material_map, Nx, Ny)
    % Initialize properties
    k = zeros(Ny, Nx); rho = zeros(Ny, Nx); c = zeros(Ny, Nx);
    T = zeros(Ny, Nx);
    % ABS
    k(material_map == 1) = 0.23;
    rho(material_map == 1) = 1040;
    c(material_map == 1) = 1.5;
    T(material_map == 1) = 21;
    % Handle
    k(material_map == 2) = 0.3;
    rho(material_map == 2) = 950;
    c(material_map == 2) = 2.0;
    T(material_map == 2) = 19;
    % Vacuum
    k(material_map == 3) = 0.0003;
    rho(material_map == 3) = 0;
    c(material_map == 3) = 0;
    T(material_map == 3) = 14;
    % Water
    k(material_map == 4) = 0.575;
    rho(material_map == 4) = 1000;
    c(material_map == 4) = 4.186;
    T(material_map == 4) = 4;
    % Ice
    k(material_map == 5) = 2.2;
    rho(material_map == 5) = 917;
    c(material_map == 5) = 2.06;
    T(material_map == 5) = -15;
    % Internal Air
    k(material_map == 6) = 0.025;
    rho(material_map == 6) = 1.2;
    c(material_map == 6) = 1.005;
    T(material_map == 6) = 8;
    % Lid
    k(material_map == 7) = 0.2;
    rho(material_map == 7) = 950;
    c(material_map == 7) = 1.75;
    T(material_map == 7) = 15;
    % Stainless Steel (Inside)
    k(material_map == 8) = 16.2;
    rho(material_map == 8) = 8000;
    c(material_map == 8) = 0.5;
    T(material_map == 8) = 5;
    % Stainless Steel (Outside)
    k(material_map == 9) = 15.2;
    rho(material_map == 9) = 8000;
    c(material_map == 9) = 0.5;
    T(material_map == 9) = 19;
    % External Air
    k(material_map == 10) = 0.025;
    rho(material_map == 10) = 1.2;
    c(material_map == 10) = 1.005;
    T(material_map == 10) = 20;
end

I have two questions. firstone is that if i change the dx,and dy to 0.5, the figure flashes with check patterns.

also i wonder why the water near inside stainless steel goes below zero. please help