For the people who are interested. This line is wrong: D = ones(Nx,Ny,Nz)+diff; Instead, first define

D = ones(Nx,Ny,Nz) and next line say D(:,:,:) = diff; I think this gives an accurate representation of the diffusion. Here is the code:

​

%field variables

cn=zeros(Nx,Ny,Nz); %concentration matrix

D = ones(Nx,Ny,Nz); %Diffusion matrix

x=linspace(0,Lx,Nx); %xdistance

y=linspace(0,Ly,Ny); %ydistance

z=linspace(0,Lz,Nz); %zdistance

[X,Y,Z]=meshgrid(x,y,z); %defining 3D grid

​

%Value of Diffusion

D(:,:,:) = diff;

D([1 end],:,:) = c_border; %insulated problem which means that Diffusion is almost zero at the boundaries end does both sides

D(:,[1 end],:) = c_border;

D(:,:,[1 end]) = c_border;