fixed bed adsorption column model-solving PDE-Langmiur isotherm

Dear all,

I am trying to get a breakthrogh curve from the equations(attached) but am not getting proper plot. I need to model a fixed bed to obtain the breakthrough curves using langmiur isotherm. I am having problem with this code but dont know where is the problem. I would appreciate any help. Thank you very much.

Cfeed = 10;     % Inlet concentration
tf = 86400;      % End time
nz = 200;       % Number of gateways
np = nz + 1;    % Number of gateposts
% Initial Concentrations
C = zeros(np,1); C(1) = Cfeed;  
Cs = zeros(np,1);
% Prepare data to pass to ode solver
c0 = [C; Cs];   
dt = tf/10000;
tspan = 0:dt:tf;
% Call ode solver
[t, c] = ode15s(@rates, tspan, c0);
% Plot results
plot(t, c(:,np)/Cfeed),grid
xlabel('time'), ylabel('Cb/Cfeed')
title('Figure 1: Exit Cb/Cfeed vs time')
% figure
% plot(t, c(:,2*np)),grid
% xlabel('time'), ylabel('Cs')
% title('Figure 2: Exit Cs vs time')
% figure
% plot(t, c(:,np)-c(:,2*np)), grid
% xlabel('time'), ylabel('Cb - Cs')
% title('Figure 3: Exit Cb-Cs vs time')
% 
% k = 2.5*10^-5; nf = 1.45;
% q = k*c(:,2*np).^(1/nf);  % Freundlcih equation
% figure
% plot(t,q), grid
% xlabel('time'), ylabel('q')
% title('Figure 4: Exit q vs time')
% Function to calculate rate of change of concentrations with time
function dcdt = rates(~, c)
    Dz = 3*10^-8;   % Diffusion coefficient
    u = 2*10^-5;    % Velocity
    e = 0.4;
    Kf = 3*10^-5;   
    ap = 2*10^-3;   % Particle diameter
    ep = 0.53;
    L = 0.1;        % Length
    nz = 200;       % Number of gateways
    np = nz + 1;    % Number of gateposts
    dz = L/nz;
    dr = ap/nz;
    rho = 400;
    Dp = 1.2*10^-10;
    b=0.03;
    qm=118;
    C = c(1:np);
    Cs = c(np+1:end);
    Cs(1) = 0;
    dCdt = zeros(np,1);
    dCsdt = zeros(np,1);
    dCsdr = zeros(np,1);
    
%     rhalf = zeros(np-1,1);
%     DCsDr = zeros(np,1);
     D2CsDr2 = zeros(np,1);
%     
%     rhalf(1:np-1)=(z(1:np-1)+z(2:np))/2;
    
    for i = 2:np-1
       dCsdr(i) = (2*dr)*(Cs(i+1)-Cs(i-1));
       D2CsDr2(i) = (1/dr^2)*(Cs(i+1)-2*Cs(i)+Cs(i-1));
      % dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i));  % Needed to change the sign from -ve to +ve here.
   end % i loop
    
   for i = 2:np-1
       dCdt(i) = (Dz/dz^2)*(C(i+1)-2*C(i)+C(i-1)) - u/(2*dz)*(C(i+1)-C(i-1)) - (1-e)*3*Kf/(e*ap)*(C(i)-Cs(i));
      % dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i));  % Needed to change the sign from -ve to +ve here.
      dCsdt(i) = 1/(1+rho*((1-ep)/ep)*(qm*b/((1+b*Cs(i))^2)))*Dp*(D2CsDr2(i)+((2/ap)*dCsdr(i)));
   end % i loop
   
   % For simplicity, set exit rates to those of the immediately previous nodes
    dCdt(np) = dCdt(np-1); 
    dCsdt(np) = dCsdt(np-1);
    
    % Combine rates for the function to return to the calling routine
    dcdt = [dCdt; dCsdt];
end % of rates function