% Example LEP-17-6a: Segregation model with Asymmetric RTD (Multiple Reactions)
%

tspan = [0 2.52]; % Range for the independent variable
y0 = [1 1 0 0 0 0 0 0 0 0 0]; % Initial values for the dependent variables
k1=1;
k2=1;
k3=1;
[t,y]=ode45(@(t,y)ODEfun(t,y,k1,k2,k3),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,500])
nexttile
plot(t,y(:,11),'LineWidth',1.5);
xlabel('time (mins)')
ylabel('$\bar{X}$','Interpreter','Latex')
leg=legend('X bar','Location','NorthWest');
axis([0 2.52 0 1])
grid on
leg.ItemTokenSize=[10,10];
title('Conversion profile')
nexttile
plot(t,y(:,4),t,y(:,5),t,y(:,6),t,y(:,9),t,y(:,10),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('Concentration')
leg=legend('C_A bar','C_B bar','C_C bar','C_D bar','C_E bar','Location','NorthWest');
leg.ItemTokenSize=[10,10];
title('Concentration profile')
axis([0 2.52 0 1])
grid on
z=size(y);
ca0 =  1;
for i=1:z(1,1)
   rc(i) = k1*y(i,1)*y(i,2);
   re(i) = k3*y(i,2)*y(i,7); 
   ra(i) = -k1*y(i,1)*y(i,2)-k2*y(i,1);
   rb(i) = -k1*y(i,1)*y(i,2)-k3*y(i,2)*y(i,7);
   rd(i) = k2*y(i,1)-k3*y(i,2)*y(i,7);
   Scd(i)=y(i,3)/(y(i,7)+0.00001);
   Sde(i)=y(i,7)/(y(i,8)+0.00001);
end
nexttile
plot(t,Scd(1:z(1,1)),t,Sde(1:z(1,1)),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('Selectivity')
leg=legend('S_{C/D}','S_{D/E}','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Selectivity profile')
axis([0 2.52 0 20])
grid on
nexttile
plot(t,-ra(1:z(1,1)),t,-rb(1:z(1,1)),t,-rc(1:z(1,1)),t,-rd(1:z(1,1)),t,-re(1:z(1,1)),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('Rate')
leg=legend('-r_A','-r_B','-r_C','-r_D','-r_E','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Rate profile')
axis([0 2.52 -3 3])
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,k1,k2,k3); 
ca= Yfuncvec(1);
cb= Yfuncvec(2);
cc= Yfuncvec(3);
cabar= Yfuncvec(4);
cbbar= Yfuncvec(5);
ccbar= Yfuncvec(6);
cd= Yfuncvec(7);
ce= Yfuncvec(8);
cdbar= Yfuncvec(9);
cebar= Yfuncvec(10);
Xbar= Yfuncvec(11);
% Explicit equations
cao =  1;
E1 = -2.104*t^4+4.167*t^3-1.596*t^2+0.353*t-0.004;
E2 = -2.104*t^4+17.037*t^3-50.247*t^2+62.964*t-27.402;
rc = k1*ca*cb;
re = k3*cb*cd; 
ra = -k1*ca*cb-k2*ca;
rb = -k1*ca*cb-k3*cb*cd;
rd = k2*ca-k3*cb*cd;
X=(cao-ca)/cao;
if t<=1.26
E = E1;
else
E=E2;
end
% Differential equations
dcadt = ra;
dcbdt = rb;
dccdt = rc;
dcabardt = ca*E;
dcbbardt = cb*E;
dccbardt = cc*E;
dcddt = rd;
dcedt = re;
dcdbardt = cd*E;
dcebardt = ce*E;
dXbardt=X*E;
dYfuncvecdt=[dcadt;dcbdt;dccdt;dcabardt;dcbbardt;dccbardt;dcddt;dcedt;dcdbardt;dcebardt;dXbardt];
end