На выходе программы получаю искривление теплового фронта - Matlab - Ответ 7562271

clc, clear all
% Дано
Lx = 10; Ly = 10; hx = 0.1; hy = 0.1; Nx = Lx/hx; Ny = Ly/hy;
lamda_fld = 1.35; lamda_sld = 2.25; rho_fld = 1250; rho_sld = 1550;
C_fld = 3325e03; C_sld = 2500e03; a_sld = lamda_sld/(rho_sld*C_sld);
a_fld = lamda_fld/(rho_fld*C_fld); T0 = -2.5; tau = 3600; t_end = 3600*24*31; Lq = 155e06; Tbf = 0;
delta = 0.1;
% Начальная температура
T = ones(Nx,Ny)*T0;
% Граничные условия
T(1,:) = 10;
T(Nx,:) = -2;
 
time = 0;
 
while time<t_end
    
    time = time + tau;
    % Учитываем нелинейность теплоемкости и теплопроводности
    for i = 1:Nx
        for j = 1:Ny
            if T(i,j) < (Tbf-delta)
                C(i,j) = C_sld;
                lamda(i,j) = lamda_sld;
            elseif T(i,j) >= (Tbf-delta) && T(i,j) <= (Tbf+delta)
                C(i,j) = 0.5*(C_fld+C_sld) + 0.5*Lq/delta;
                lamda(i,j) = lamda_sld + 0.5*(lamda_fld-lamda_sld)/delta*(T(i,j) - (Tbf - delta));
            else
                C(i,j) = C_fld;
                lamda(i,j) = lamda_fld;
            end
        end
    end
    % Здаем нулевой теплопоток на границах
    T(:,1) = T(:,2);
    T(:,Ny) = T(:,Ny-1);
    % Считаем уравнение теплопроводности
    for i = 2:Nx-1
        for j = 2:Ny-1
            T(i,j) = T(i,j) + tau/hx^2*(lamda(i,j)/C(i,j))*(T(i+1,j) - 2*T(i,j) + T(i-1,j))+...
                tau/hy^2*(lamda(i,j)/C(i,j))*(T(i,j+1) - 2*T(i,j) + T(i,j-1));
        end
    end
end
X = 0:Lx/(Nx-1):Lx; Y = 0:Ly/(Ny-1):Ly;
contourf(X, Y, T) % визуализируем
% colorbar, 
grid on % закрепляем сетку и цвет.шкалу
title('Температурное поле массива грунта на 30 год эксплуатации (на 15 сентября)') %название
xlabel('Ширина, м'); ylabel('Глубина, м'); zlabel('Температура, deg[C]') %название осей
axis auto, caxis([-5 5]) %X и Y оси автоматически, T - в заданных пределах