% LEP-10 -Algorithm for decaying catalyst in a moving bed reactor
%

Wspan = [0 500]; % Range for the independent variable W
y0 = [0;1;1]; % Initial values for the dependent variables i.e X,a,p
k =0.1;
KA =1.5;
kd =0.75;
FA0 =10;
PA0= 20;
Us=2;
alpha=0.0019;
[w,y]=ode45(@(w,y)ODEfun(w,y,k,KA,kd,FA0,PA0,Us,alpha),Wspan,y0);
tiledlayout(1,2)
set (gcf,'Position',[0,0.1,800,300])
nexttile
plot(w,y(:,1),w,y(:,2),'LineWidth',2)
ylabel('X, p');
xlabel('W(kg)');
leg=legend('X','p');
leg.ItemTokenSize=[10,10];
axis([0 500 0 1])
title('Conversion and pressure drop profiles')
grid on
nexttile
plot(w,y(:,3),'LineWidth',2)
ylabel('Activity');
xlabel('W(kg)');
leg=legend('a');
leg.ItemTokenSize=[10,10];
axis([0 500 0 1])
title('Activity profile')
grid on

%

function dYfuncvecdZ = ODEfun(w,Y,k,KA,kd,FA0,PA0,Us,alpha); 
X=Y(1);
p=Y(2);
a=Y(3);
% Explicit equations
PA= PA0*(1 - X)*p;
PB= PA0*X*p;
rA= -a*k*PA/(1 + KA*PA);
% Differential equations
dXdw = -rA/FA0;
dadw=-kd*(a^2)*PB/Us;
dpdw=-alpha/2*p;
dYfuncvecdZ = [dXdw;dpdw;dadw]; 
end
