Can we solve 3 pde simultaneously using pdepe with respect to space and time??

I have a system of 3 equations with 3 variables and corresponding boundary conditions. I am using pdepe to solve these. Is it possible? Can anyone help me correcting the code I have written? Sharing the equation in an image attached. 3rd equation is : dQp/dt=Keq*Cp*(lambda-(sigma+charge)*Qp)^charge-Qp*Csin^charge

%% Main function
function sma_porediffusion
%% globalising variables and constants
global v x F Ac Dcoe ap Cin charge lambda epsilon k_eq sigma t m L dc N Run_time CV Csin;
m=0;

%% specifying system parameters
ap=.0023; %particle size
epsilon_e= 0.7; %porosity
epsilon_p=0.8; %
L =5.2; %cm
dc= 0.7; %cm
F=0.5;%ml/min
Ac=3.142*(dc/2)^2; %cm2
v=F/Ac;%cm/min
Dcoe=4.43*10^-3;%cm^2/min
charge=2.14;
lambda= 1833; % capacity factor in mM
Cin=0.0085; %mg/ml
k_eq= 0.12; % equilibrium coefficient = kads/kdes
sigma= 49; % steric hinderence factor
N=50; % meshing intervals
Run_time =1; % min
Csin= 1; %mM mobile phase salt concentration
% Csout= 400;%mM end salt concentration
% dq/dt=keq*C*(lambda-(sigma+charge)*q)^charge-q*(Csalt^charge)
% Area= pi*(dc^2)/4; % area of cross section cm^2

%%
x=linspace(0,L,N);
t=linspace(0,Run_time,N);
sol=pdepe(m,@sma_pdecolumnpe,@sma_pdecolumnic,@sma_pdecolumnbc,x,t);

%% Extract the first solution component as conc.
conc1 = sol(:,:,1);
conc2=conc1./Cin;
%surf(x,t,conc2)

%% A solution profile for conc. vs length.
% figure, plot(x,conc1(end,:))
% title(strcat('Solution at t = ', num2str(t)))
% xlabel('Distance x')
% ylabel('Concentration (C)')
% 
% %Plot conc. vs. time
% figure, plot(t,conc1)
% title('conc vs time: breakthrough conc1')
% xlabel('Time (min)')
% ylabel('Concentration (mg/ml)')

figure, plot(t,conc1(:,end))
title('conc vs time: breakthrough')
xlabel('Time (min)')
ylabel('Concentration (mg/ml)')

%% Inlet conditions
function u0 = sma_pdecolumnic(x)
if x==0
    u0=[Cin;0;0];
else
    u0=[0;0;0];
end
end
%% describing equations
function [c,f,s] = sma_pdecolumnpe(x,t,u,DuDx)
%Defining variables
%global v Cin charge lambda epsilon k_eq sigma k_a k_d Csin;
c=[1;1;1];
f=[Dcoe;0;0].*DuDx;
%dq/dt=adsorption coefficient*c(ionic capacity-(steric
%factor+characteristic charge)*q)^ char charge - desorption coeff*q*Csaltin^
%char charge

sq=k_eq.*u(1).*(lambda-(sigma+charge).*u(2))^charge- u(2)*Csin^charge;
%disp(sq);
s=[-v*DuDx(1)-(1-epsilon_e)*k_eq*ap*(u(1)-u(2));(1-epsilon_p)*sq+k_eq*ap*(u(1)-u(2));sq];
end
%% boundary conditions
function [pl,ql,pr,qr] = sma_pdecolumnbc(xl,ul,xr,ur,t)
%global Cin;
pl = [ul(1)-Cin;0;0];
ql = [0;0;0]; 
pr = [0;0;0];
qr = [1;1;1];
end
end

Thanks for your inputs