% LEP 10-5 : Catalytic Cracking in a Moving Bed Reactor
%

Wspan = [0 22]; % Range for the independent variable W
y0 = [1,0]; % Initial values for the dependent variables i.e a and X
Us=10;
kd=0.72;
Fao=30;
Cao=0.075;
k=600;
[w,y]=ode45(@(w,y)ODEfun(w,y,Us,kd,Fao,Cao,k),Wspan,y0);
tiledlayout(1,2)
set (gcf,'Position',[0,0.1,800,300])
nexttile
plot(w,y(:,1),w,y(:,2),'LineWidth',1.5)
ylabel('a, X');
xlabel('W (kg)');
leg=legend('a','X');
leg.ItemTokenSize=[10,10];
axis([0 22 0 1])
title('Conversion and Activity profiles')
grid on
z=size(y);
for i=1:z(1,1)
    Ca(i)=Cao*(1-y(i,2));
    raprime(i)=-k*Ca(i)^2;
end
nexttile
plot(w,-raprime(1:z(1,1)),'LineWidth',1.5)
ylabel('Rate (mol/kg.cat.min)');
xlabel('W (kg)');
leg=legend('-r_A^\prime');
leg.ItemTokenSize=[10,10];
axis([0 22 0 5])
title('Rate profile')
grid on

%

function dYfuncvecdW = ODEfun(W,Y,Us,kd,Fao,Cao,k); 
a = Y(1);
X=Y(2);
% Explicit equations
Ca=Cao*(1-X);
raprime=-k*Ca^2;
% Differential equations
dadW= -kd*a/Us; 
dXdW=a*(-raprime)/Fao;
dYfuncvecdW = [dadW;dXdW]; 
end
