% Example LEP-WEB17-1h: Maximum Mixedness Model with Asymmetric RTD (Multiple Reactions)
%

zspan = [0 2.52]; % Range for the independent variable
y0 = [1 1 0 0 0 0.99]; % Initial values for the dependent variables
k1=1;
k2=1;
k3=1;
[z,y]=ode45(@(z,y)ODEfun(z,y,k1,k2,k3),zspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,500])
nexttile
plot(z,y(:,1),z,y(:,2),z,y(:,3),z,y(:,4),z,y(:,5),'LineWidth',1.5);
xlabel('z (mins)')
ylabel('Concentration')
leg=legend('C_Abar','C_Bbar','C_Cbar','C_Dbar','C_Ebar','Location','NorthEast');
axis([0 2.52 0 1])
grid on
leg.ItemTokenSize=[10,10];
title('Concentration profile')
p=size(y);
cao =  1;
for i=1:p(1,1)
   X(i)=1-y(i,1)/cao;  
   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,4); 
   rc(i) = k1*y(i,1)*y(i,2); 
   rd(i) = k2*y(i,1)-k3*y(i,2)*y(i,4); 
   re(i) = k3*y(i,2)*y(i,4);
   Scd(i)=y(i,3)/(y(i,4)+0.00001);
   Sde(i)=y(i,4)/(y(i,5)+0.00001);
end
nexttile
plot(z,X(1:p(1,1)),'LineWidth',1.5)
xlabel('z(mins)')
ylabel('Conversion')
leg=legend('Xbar','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Conversion profile')
axis([0 2.52 0 1])
grid on
nexttile
plot(z,-ra(1:p(1,1)),z,-rb(1:p(1,1)),z,-rc(1:p(1,1)),z,-rd(1:p(1,1)),z,-re(1:p(1,1)),'LineWidth',1.5)
xlabel('z(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
nexttile
plot(z,Scd(1:p(1,1)),z,Sde(1:p(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

%

function dYfuncvecdz = ODEfun(z,Yfuncvec,k1,k2,k3); 
ca= Yfuncvec(1);
cb= Yfuncvec(2);
cc= Yfuncvec(3);
cd= Yfuncvec(4);
ce= Yfuncvec(5);
F= Yfuncvec(6);
% Explicit equations
ra = -k1*ca*cb-k2*ca;
rb = -k1*ca*cb-k3*cb*cd; 
rc = k1*ca*cb; 
rd = k2*ca-k3*cb*cd; 
re = k3*cb*cd; 
cao = 1;
cbo = 1;
cco = 0;
cdo = 0;
ceo = 0 ;
lam = 2.52 - z;
E1=-2.104*(2.52 - z).^4 + 4.167*(2.52 - z).^3 - 1.596*(2.52 - z).^2 + 0.353*(2.52 - z) - 0.004;
E2=-2.104*(2.52 - z).^4 + 17.037*(2.52 - z).^3 - 50.247*(2.52 - z).^2 + 62.964*(2.52 - z) - 27.402;
if lam<=1.26
    E=E1;
else
  E=E2;
end 
EF=E/(1-F);
% Differential equations
dcadz =   -(-ra+(ca-cao)*EF) ;
dcbdz =   -(-rb+(cb-cbo)*EF) ; 
dccdz =  -(-rc+(cc-cco)*EF) ; 
dcddz =  -(-rd+(cd-cdo)*EF) ; 
dcedz =  -(-re+(ce-ceo)*EF) ; 
dFdz =  -E;
dYfuncvecdz=[dcadz;dcbdz;dccdz;dcddz;dcedz;dFdz];
end