% Defining grid in the horizontal direction
xmin=0.0; xmax=3000.0; xnum=99;
xstp=(xmax-xmin)/(xnum-1);
ymin=0.0; ymax=3000.0; ynum=99;
ystp=(ymax-ymin)/(ynum-1);

k=3.0; % Thermal conductivity, W/m/K
ro=3000.0; %Density, kg/m^3
cp=1000.0; % Isobaric heat capacity, J/kg
kappa=k/ro/cp; % Thermal diffusivity, m^2/s

timestep=0.75*((xstp^2+ystp^2)/2.0)/kappa/3.0; % Crucial parameter: Explicit timestep
% Defining Initial Temperature structure
tmin=500;
tmax=1000;
for xn = 1:xnum
for yn = 1:ynum
T0(xn,yn)=tmin; % Background Temperature, K
if ((xn-1)/(xnum-1)>0.25 & (xn-1)/(xnum-1)<0.75 & (yn-1)/(ynum-1)>0.25 & (yn-1)/(ynum-1)<0.75)
T0(xn,yn)=tmax; % Rectangular Hot body in the middle
end
end
end

% Solving temperature equation by explicit method
for cycle = 0:25
for xn = 1:xnum
for yn = 1:ynum
% Boundary nodes:
if (xn==1 | xn==xnum | yn==1 | yn==ynum)
T1(xn,yn)=tmin; % T=const
else
kappaX=kappa;
kappaY=kappa;
T1(xn,yn)=T0(xn,yn)+timestep*kappaX*( T0(xn-1,yn)-2.0.*T0(xn,yn)+T0(xn+1,yn) )/(xstp*xstp); %x contribution
T1(xn,yn)=T1(xn,yn)+timestep*kappaY*( T0(xn,yn-1)-2.0.*T0(xn,yn)+T0(xn,yn+1) )/(ystp*ystp); %y contribution
end
end
end
T0=T1;% Reloading solutions to T0
% Creating vectors for x and y axes
x = xmin:xstp:xmax;
y = ymin:ystp:ymax;
x=x/1000.;% renormalizes m in km, only for plotting
y=y/1000.;% renormalizes m in km, only for plotting
% Visualizing results
% Ploting RO as colormap
figure(1);

surf(x,y,T0);
view (0,90);
caxis([tmin tmax]);
colorbar;
title(['Temperature, K For time step ',num2str(cycle)])
xlabel('X, km')
ylabel('Y, km')
end