Multispecies Diffusion - Fick's 2nd Law (PDEs)

I am trying to simulate an unsteady state multispecies diffusion scenario, in Matlab using a system of PDEs (Fick's Second Law) with discontinuous initial conditions. The geometry at t=0 is shown below.

Initially (at t=0), in region one only species 1 is present, and in region 2 only species 2, and 3 are present. As t increases, all species begin to diffuse, until steady state is reached. The boundary conditions I am using, are insulated boundary conditions, but I seem to be missing a crucial element that couples the PDEs, maybe? Not sure.

The system of PDEs I am attempting to solve are the simplest versions of Fick's second law and are as follows:

When I attempt to solve the system of PDEs using pdepe built-in function in MATLAB, I keep getting zero concentrations for all species throughout the domain at all time points. I am unsure how to couple the PDEs. Should there be another term in the system of PDEs? What am I missing?

L = 15;
x = linspace(0,L,30);
t = [linspace(0,3600,30)];
m = 0;
pdepe(m,@pdefun,@icfun,@bcfun,x,t)

function [c,f,s] = pdefun(x,t,C,dCdx)
D = [1.79e-14;2.51e-14;1.84e-14];
c = [1;1;1];
f = D.*[dCdx(1);dCdx(2);dCdx(3)];
s = [0;0;0];
end

function C0 = icfun(x,t)
    if x<=9  
        C0 = [1056;0;0];
    else
        C0 = [0;254;1343];
    end 
end

function [pL,qL,pR,qR] = bcfun(xL,CL,xR,CR,t)
pR = [0;2.51e-14;1.84e-14];
qR = [1;1;1];
pL= [1.79e-14;0;0];
qL = [1;1;1];
%end
end

Answer 1

[Torsten]

Check the consistency of the length and time units for your constants L, t, D and C0.

L = 15;
x = linspace(0,L,3000);
%t = linspace(0,3600,30);
t = linspace(0,360000,30);
m = 0;
sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);
u = sol(:,:,1);
plot(x,[u(10,:);u(20,:);u(30,:)])

function [c,f,s] = pdefun(x,t,C,dCdx)
%D = [1.79e-14;2.51e-14;1.84e-14];
D = 1000000000*[1.79e-14;2.51e-14;1.84e-14];
c = [1;1;1];
f = D.*[dCdx(1);dCdx(2);dCdx(3)];
s = [0;0;0];
end

function C0 = icfun(x,t)
    if x<=9  
        C0 = [1056;0;0];
    else
        C0 = [0;254;1343];
    end 
end

function [pL,qL,pR,qR] = bcfun(xL,CL,xR,CR,t)
%pR = [0;2.51e-14;1.84e-14];
pR = [0;0;0];
qR = [1;1;1];
%pL= [1.79e-14;0;0];
pL = [0;0;0];
qL = [1;1;1];
%end
end

[Rawatee Maharaj-Sharma]

Makes sense. Thank you, seems like I overlooked difference in units!