problem at solving coupled PDE-ODE with reaction between 4 species

I have solved a time dependent coupled PDE-ODE system for 4 compound where there are some rection between them. my probem is that when i add reaction terms in my model, i got wrong result.the first coumpound should decrease and three others should increase but first one remains fixed and others decrease. that is very strange for me. i am sure i wrote correct expression for reactions but the result is reverse. this is my code:

     clc
  clear 
  format long
  tf=60*1*60;
  nt=100;
  dt=tf/nt;
  tspan = 0:dt:tf;
  nx = 50;
  %%%%%%%parameters%%%%
  Dab = 11.3*(1e-6);%
  ux=0.5;
  L=2e-3*(5/95);
  Cin=[ 286.4879091 3.508114893  3.769381628  9.247436195]*1e-3*2.95;
  yout=[233.9457273  5.289110609  26.08221981  31.552135]*1e-3*2.95;
  K=[11.55507758  1.841304056  561.4994587  60.69688482  1.959738497];
  Ncmp=4;
  for i=1:8*nx
      y0(1,i)=0; 
  end

options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8, 'InitialStep', 0.0001);
[t,y] = ode15s(@fun,tspan,y0,options,nx,Ncmp,L,Dab,Cin,ux,K);
[s1,s2] = size(t);
plot(t/60,y(1:s1,2*nx-1))
hold on
plot(t/60,y(1:s1,4*nx-1))
hold on
plot(t/60,y(1:s1,6*nx-1))
hold on
plot(t/60,y(1:s1,8*nx-1))
hold on
plot(25,yout,'r*')
function dydt=fun(~,y,nx,Ncmp,L,Dab,Cin,ux,K)
kg=0.66;      
area=4; 
eps=0.65;
rho=4.23*1e6;

B=rho;
E=1;  %<---E can be changed

dx=L/(nx+1);
dydt=zeros(8*nx,1);
k1=K(1);
k1a=K(2);
k2=K(3);
k3=K(4);
k4=K(5);
for j=1:Ncmp
   for i=1:2*nx
   C(j,i)=y(i+2*nx*(j-1));     % concentration of j in node i, 
   end
end
%%%%%%%%%%%Reactions%%%%%%%%
for i=2:2:2*nx
    rd1(i)=k1*k1a*C(1,i)/(1+k1a*C(1,i));
    rd2(i)=k2*C(2,i);
    rd3(i)=k3*C(3,i);
    rd4(i)=k4*C(4,i);
    react(1,i)=-rd1(i);
    react(2,i)=rd1(i)-rd2(i);
    react(3,i)=rd2(i)-rd3(i);
    react(4,i)=rd3(i)-rd4(i);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
end
for j=1:Ncmp
i=1;
        C0(j)=Cin(j);
        dcdx(j,i)=(C(j,i+2)-C0(j))/2/dx;
        d2cdx2(j,i)=(C(j,i+2)-2*C(j,i)+C0(j))/dx/dx;
        dcdT(j,i)= Dab*d2cdx2(j,i)-ux*dcdx(j,i)-kg*area*((1-eps)/eps)*(B*C(j,i)-C(j,i+1));
        dcdT(j,i+1)=kg*area*(B*C(j,i)-C(j,i+1))+E*react(j,i+1);    %<---E can be changed 

for i = 3:2:2*nx-3
        dcdx(j,i)=(C(j,i+2)-C(j,i-2))/2/dx;
        d2cdx2(j,i)=(C(j,i+2)-2*C(j,i)+C(j,i-2))/dx/dx;
        dcdT(j,i)= Dab*d2cdx2(j,i)-ux*dcdx(j,i)-kg*area*((1-eps)/eps)*(B*C(j,i)-C(j,i+1));
        dcdT(j,i+1)=kg*area*(B*C(j,i)-C(j,i+1))+E*react(j,i+1);          %<---E can be changed
end
i = 2*nx-1;
C(j,i+2)=C(j,i)-(C(j,i-2)-C(j,i))/3;
         dcdx(j,i)=(C(j,i+2)-C(j,i-2))/2/dx;
        d2cdx2(j,i)=(C(j,i+2)-2*C(j,i)+C(j,i-2))/dx/dx;
        dcdT(j,i)= Dab*d2cdx2(j,i)-ux*dcdx(j,i)-kg*area*((1-eps)/eps)*(B*C(j,i)-C(j,i+1));
        dcdT(j,i+1)=kg*area*(B*C(j,i)-C(j,i+1))+E*react(j,i+1);         %<---E can be changed
end
 for j=1:Ncmp
    for i = 1:2*nx
    dydt(i+2*nx*(j-1)) = dcdT(j,i);
    end
 end      

end

[Moji]

when i use this code without reaction terms ( * react(j,i)*)the curve and result is reasonable(fig 1) but with applying reaction terms(fig 2) it goes wrong. i can not find where is my mistake. i would be grateful if anybody could help me.

[Torsten]

Please plot concentrations over length for several times instead of concentrations over time for some fixed position. It makes it much easier to interpret the results.

[Moji]

@Torsten thank you for your reply. i put the figure for different time (t=2,100) over the length(at t=1 all concentration are zero as my initial caondition). as you can see concentration of othe compounds(B,C,D) doesnt change at all. it means other compounds doesn't create but based on my reaction expression they should be created and their concentrations should be incresed.

[Moji]

if i remve the reaction terms the result would be like this figure. it is wierd, i dont know why concentrations of B,C,D decreases instead of increase, when reaction terms are added.

[Torsten]

The usual consequence would be to output react(j,i).

[Moji]

i couldn't get your point, could you please explain a little more?

[Torsten]

You should check whether react(j,i) is negative for j=2,3,4 (because you see decreasing concentrations for B,C and D). And if yes: whether the K-array is reasonable.

[Moji]

i change it for react(j,i), it works now. i really really really really thank you. but the reaction is mentioned in ODE (solid phase), why this happened?

[Torsten]

i change it for react(j,i), it works now.

I don't understand what you did.

[Moji]

in this line :

dcdT(j,i+1)=kg*area*(B*C(j,i)-C(j,i+1))+E*react(j,i+1); 
i changed react(j,i+1) to react(j,i). and also i changed the for loop in below lines. i used for i=1:2:2*nx instead of for i=2:2:2*nx:

for i=2:2:2*nx
    rd1(i)=k1*k1a*C(1,i)/(1+k1a*C(1,i));
    rd2(i)=k2*C(2,i);
    rd3(i)=k3*C(3,i);
    rd4(i)=k4*C(4,i);
    react(1,i)=-rd1(i);
    react(2,i)=rd1(i)-rd2(i);
    react(3,i)=rd2(i)-rd3(i);
    react(4,i)=rd3(i)-rd4(i);