Уравнение теплопроводности

close all
clear all
clc
 
dx=1; dy=1; dt=.1; 
x=0:dx:199; y=0:dy:199; lx=length(x); ly=length(y);
 
u=zeros(lx,ly);
un=zeros(lx,ly);
u(1:30,1:end)=-5; %temperature of depth 1
u(30:75,1:end)=-4; %temperature of depth 2
u(75:150,1:end)=-2; %temperature of depth 3
u(155:end,1:end)=0; %temperature of depth 4
u(1:end,1)=24; %temperature left
 
for tt=0:500; %timestep    
    for jj=2:ly-1;
        for ii=2:lx-1;
            if ii && ii<=30 %layer 1
                k=1.86;
                cp=1.8855;
                rho=1.400;
            elseif ii>30 && ii<=55 %layer 2
                k=1.34;
                cp=2.304;
                rho=0.500;
            elseif ii>125 && ii<=55 %layer 3               
                k=0.81;
                cp=3.210;
                rho=0.500;
            else 
                k=15.5; %layer 4
                cp=0.450;
                rho=8.900;              
            end
            if tt && tt<=25 
                u(200,1:end)=17; %temperature of time 1
            elseif tt>26 && tt<=89
                u(200,1:end)=-15; %temperature of time 2
            elseif tt>90 && tt<=150               
                u(200,1:end)=25; %temperature of time 3
            else 
                u(200,1:end)=-5; %temperature of time 4          
            end
            u(ii,jj)=u(ii,jj)+(k/(cp*rho))*(dt/dx^2)*...
            (u(ii+1,jj)-2*u(ii,jj)+u(ii-1,jj))+...
            (k/(cp*rho))*(dt/dy^2)*(u(ii,jj+1)+...
            u(ii,jj-1)-2*u(ii,jj));  %heat transfer equation
        end
        un=u;
    end
    %imagenation
    imagesc(un)
    set(gca,'ydir','normal')
    colorbar;
    caxis([-10 10])
    pause(eps)
    drawnow;
    hold on
end