Boundary conditions in finite element problem

Good afternoon everyone,

I am trying to use the finite element analysis of the PDE toolbox for a simple diffusion problem. The geometry consists of a sphere within a sphere, and they are concentric. The inner sphere has a certain substance (Initial concentration = 1) , but the outer doesn't (I.C = 0). Furthermore, both domains have different diffusion coefficients (inner is faster than outer). The boundary condition on the exterior of the outer sphere is set to 0. So the substance is supposed to diffuse from the inner to the outer. Both domains obey Fick's law, the difference is the diffusion coefficient within each domain.

So, as far as I understand, I can do this 3D shape as a 2D axisymmetric problem. I have worked for a while in my code, but it is not working properly. For instance: at t=0, I see in my plot that the basically both domains have a concentration of 1 (instead of the inner being 1 and the outer being 0); and also for certain time points I am getting negative values, which is not expected. Any comments/suggestions?

%% Step 1. Geometry set up
myPDE = createpde(2); %I believe is 2, since I have 2 PDE equations, 1 for each domain
geom1 = decsg([1 0 0 0.00000635000000000000; 1 0 0 0.00000510000000000000]'); %Concentric circles
geometryFromEdges(myPDE,geom1);
figure (1)
pdegplot(myPDE,'FaceLabels','on','EdgeLabels','on')
axis equal


%% Step 2. Specify PDE coefficients
%Following three coefficients are interpreted as diagonal matrix if they are scalars
d=[1;1];
m=0;
a=0;
d1=1e-14; %Diffusion coefficient in domain 1
d2=1e-15; %Diffusion coefficient in domain 2
c=[d1;d2];
f=[0; 0];
specifyCoefficients(myPDE,'m',m,'d',d,'c',c,'a',a,'f',f); %Parabolic

%% Step 3. B.C
applyBoundaryCondition(myPDE,'dirichlet','edge',[1,2,3,4,5],'u',[0,0]); %Concentration in the outermost part is 0

%% Step 4. Set I.C
c0 = [1;0]; %Initial concentration in the inner domain; I.C in the outer domain
setInitialConditions(myPDE,c0);

%% Step 5. Create mesh
generateMesh(myPDE);
figure (2)
pdeplot(myPDE)

%% Step 6. Solve
timelist = 0:1000:170000;
result = solvepde(myPDE,timelist); 

u = result.NodalSolution;
pdeplot(myPDE,"XYDATA",u(:,1));

[Torsten]

Why don't you use pdepe on the problem

du/dt = 1/r^2 * d/dr(r^2*D(r)*du/dr)

?