Problem with time plot from partial derivative equation solved by method of lines

Hello, I have a problem ,I tried to solve an pde equation using method of lines but appears a warning when time increases that dont let me raise the time and has a limit at time 7.5 h , dont know to solve it . I show the equations to better see the problem.Thanks a lot.

desorcion_final(0.4,0.5,2,0.02,20,1520,8,50.63,3.76,100,"si")

Warning: Failure at t=2.777457e+04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.820766e-11) at time t.

Warning: Imaginary parts of complex X and/or Y arguments ignored.

ans = 
   0.7441 - 0.0000i   7.5932 + 0.0000i

%tf means the time in hours
function [a]=desorcion_final(L,eps,u,k_f,c0,rhop,tf,Kf,n,q_inicial,sn)

Nt=60;
tf=tf*3600;
t=linspace(0,tf,Nt);


Nz=70;
z=linspace(0,L,Nz);
dz=z(2)-z(1);


q0star=Kf*c0^(1/n);

%Initial conditions
ICA=ones(1,Nz)*c0;
ICB=ones(1,Nz)*q_inicial;
IC=[ICA ICB];


[t, y]=ode15s(@fun_pde,t,IC,[],Nz,eps,n,Kf,k_f,u,rhop,dz);


%define value
cc=y(:,1:Nz);
qq=y(:,Nz+1:2*Nz);

%recalculate new limit conditions
cc(:,1)=0;
cc(:,end)=cc(:,end-1); %dC/dz=0

%plotting
 
cp=cc(:,end)./c0;
qp=qq(:,end-1)/q_inicial;

if sn=="si"
plot(t/3600,(1-qp))
title('Curve desorption conc')
xlabel('time(hora)')
ylabel('Convertion')
end
%% 
%Return 

a(1)=(1-qp(end));%conversion;
a(2)=t(end)/3600; %tiempo;

%% 

% PDE function

function dydt=fun_pde(t,y,Nz,eps,n,KF,k_f,u,rhop,dz)

dcdt=zeros(Nz,1);
dqdt=zeros(Nz,1);

c=y(1:Nz);
q=y(Nz+1:2*Nz);
%BC
c(1)=0; %BC1
c(end)=c(end-1);%BC2
%interior

for i=2:Nz-1
    qstar(i)=KF.*c(i).^(1/n);%Freundlich isotherm
    dqdt(i)=k_f.*(qstar(i)-q(i));
    dcdz(i)=(c(i+1)-c(i-1))./2./dz;
    dcdt(i)=-u*dcdz(i)-rhop*((1-eps)./eps).*dqdt(i);
  
end
dydt=[dcdt;dqdt];
end
end

Answer 1

[Torsten]

It's "Conversion", not "Convertion".

And setting a boundary condition at z = L is incorrect because you only have convection, no diffusion.

desorcion_final(0.4,0.5,2,0.02,20,1520,8,50.63,3.76,100,"si")

ans = 1×2

0.7497 8.0000

%tf means the time in hours
function [a]=desorcion_final(L,eps,u,k_f,c0,rhop,tf,Kf,n,q_inicial,sn)

Nt=60;
tf=tf*3600;
t=linspace(0,tf,Nt);


Nz=70;
z=linspace(0,L,Nz);
dz=z(2)-z(1);


q0star=Kf*c0^(1/n);

%Initial conditions
ICA=ones(1,Nz)*c0;
ICB=ones(1,Nz)*q_inicial;
IC=[ICA ICB];


[t, y]=ode15s(@fun_pde,t,IC,[],Nz,eps,n,Kf,k_f,u,rhop,dz);


%define value
cc=y(:,1:Nz);
qq=y(:,Nz+1:2*Nz);

%recalculate new limit conditions
cc(:,1)=0;
%cc(:,end)=cc(:,end-1); %dC/dz=0

%plotting
 
cp=cc(:,end)./c0;
qp=qq(:,end-1)/q_inicial;

if sn=="si"
plot(t/3600,(1-qp))
title('Curve desorption conc')
xlabel('time(hora)')
ylabel('Convertion')
end
%% 
%Return 

a(1)=(1-qp(end));%conversion;
a(2)=t(end)/3600; %tiempo;

%% 

% PDE function

function dydt=fun_pde(t,y,Nz,eps,n,KF,k_f,u,rhop,dz)

dcdt=zeros(Nz,1);
dqdt=zeros(Nz,1);

c=y(1:Nz);
q=y(Nz+1:2*Nz);
%BC
c(1)=0; %BC1
%c(end)=c(end-1);%BC2
%interior

for i=2:Nz
    qstar(i)=KF.*c(i).^(1/n);%Freundlich isotherm
    dqdt(i)=k_f.*(qstar(i)-q(i));
    if i < Nz
      dcdz(i)=(c(i+1)-c(i-1))/2/dz;
    else
      dcdz(i)=(c(i)-c(i-1))/dz; 
    end
    dcdt(i)=-u*dcdz(i)-rhop*((1-eps)./eps).*dqdt(i);
  
end
dydt=[dcdt;dqdt];
end
end