Fixed-bed adsorption - PDEPE?

Dear support team!

I am trying so solve a convection dispersion equation with pdepe and I am only getting constant outlet concentration which is not correct. I am having problem figuring it out where I am wrong?

I think the equation should get solved by pdepe according to matlab (https://se.mathworks.com/help/matlab/ref/pdepe.html). But i am not sure if I am correct in defining the "c" parameter?

Deeply appreciate any help.

The equations are:

And here is the code:

function DiffusionConvectionReaction
% Solving advection-reaction-dispersion equation for a packed bed reactor with
% pdepe matlab function
clc
clear
close all

global u D C0 x t rho eps
% define the reactor parameters
L = 0.2;                                            % reactor length (m)
Pe = 10;                                            % Peclet number
u = 0.06;                                           % average upflow velocity (m/s)
D = u*L/Pe;                                         % axial dispersion (m2/s)
C0 = 1e-1;                                          % initial concentration (kg/m3)
rho = 840;                                          % density (kg/m3)
eps = 0.3;                                         % porosity

% solution of a single pde
m = 0;                                              % assuming slab symmetry (-)
x = linspace(0,L,10);                               % reactor's length
t = linspace(0,20,100);                             % operation time (s)

sol = pdepe(m,@pdefun,@icfun,@bcfun,x,t);
c = sol;

figure;
plot(t,c(:,end));

 %-------------------------------------------------------------------------
 % define the flux and source term
    function [g,f,s] = pdefun(x,t,c,DcDx)
    qmax = 4.3e-3;
    b = 0.84e3;
    g = 1+ rho*qmax*b/(eps*(1+b*c));
    f = D*DcDx;
    s = -u*DcDx;
    end
  %------------------------------------------------------------------------
  % define initial condition
    function c0 = icfun(x)
        
    c0 = C0;                           
    end
 %-------------------------------------------------------------------------
 % define boundary conditions
    function [pl,ql,pr,qr] = bcfun(xl,cl,xr,cr,t)
        
    pl = u*(C0-cl);
    ql = 1;
    pr = 0;
    qr = 1/D;
    end
end

[Bill Greene]

Actually, your solution equals the inital value at all points for all time. This is expected from your equation and boundary conditions with a constant initial condition.

[Hasti]

Thank you very much for your reply. Yes, If the IC is changed to zero, as given in the equations, then the outlet concentration is not constant and is increasing with time (and this is what I am expecting from the model).