Matlab Code

%% code
X=52; %1=2mm
Y=52;
Z=21;

T=60*12; %24h
%clear value of u
u=zeros(T,X,Y,Z);
%diffusion coefficient
D=0.027*4/5;

for i=6:1:46
for j=6:1:46

u(1,i,j,1)=8*182.5e-02;
end
end


for i=1:X;
for k=1:Z;
u(1,i,1,k)=0;
u(1,i,Y,k)=0;
end
end
for j=1:Y;
for k=1:Z;
u(1,1,j,k)=0;
u(1,X,j,k)=0;
end
end

for t=2:1:T
for i=1:X;
for j=1:Y;
for k=1:Z;
u(t,1,j,k)=0;
u(t,i,1,k)=0;
u(t,X,j,k)=0;
u(t,i,Y,k)=0;
end
end
end

for i=2:1:X-1;
for j=2:1:Y-1;
for k=1:Z;
if k==1;
u(t,i,j,k)=u(t-1,i,j,k)+D*(u(t-1,i+1,j,k)-2*u(t-1,i,j,k)+u(t-1,i-1,j,k)+u(t-1,i,j+1,k)-2*u(t-1,i,j,k)+u(t-1,i,j-1,k)+u(t-1,i,j,k+1)-2*u(t-1,i,j,k)+0);
else
if k==Z;
u(t,i,j,k)=u(t-1,i,j,k)+D*(u(t-1,i+1,j,k)-2*u(t-1,i,j,k)+u(t-1,i-1,j,k)+u(t-1,i,j+1,k)-2*u(t-1,i,j,k)+u(t-1,i,j-1,k)+0-2*u(t-1,i,j,k)+u(t-1,i,j,k-1));
else
u(t,i,j,k)=u(t-1,i,j,k)+D*(u(t-1,i+1,j,k)-2*u(t-1,i,j,k)+u(t-1,i-1,j,k)+u(t-1,i,j+1,k)-2*u(t-1,i,j,k)+u(t-1,i,j-1,k)+u(t-1,i,j,k+1)-2*u(t-1,i,j,k)+u(t-1,i,j,k-1));
end
end
end
end
end
end

figure

w=zeros(X,Y);
v=zeros(X,Y);
for x=1:X;
for y=1:Y;
w(x,y)=u(T,x,y,1);

end
end
figure(1)
pcolor(w);


for z=1:Z;
for y=1:Y;
v(y,z)=u(T,26,y,z);

end
end
figure(2)
pcolor(v)

for x=1:X;
for y=1:Y;
ww(x,y)=u(T,x,y,Z);
end
end

figure(3)
pcolor(ww)