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")
%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
