% Example LEP-WEB17-1g: Segregation Model with Bimodal distribution (Multiple Reactions)
%

tspan = [0 6]; % 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(:,4),t,y(:,5),t,y(:,6),t,y(:,9),t,y(:,10),'LineWidth',1.5);
xlabel('time (mins)')
ylabel('Concentration')
leg=legend('C_Abar','C_Bbar','C_Cbar','C_dbar','C_ebar','Location','NorthEast');
axis([0 6 0 1])
grid on
leg.ItemTokenSize=[10,10];
title('Concentration profile')
z=size(y);
cao =  1;
for i=1:z(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,7);  
   rc(i) =  k1*y(i,1)*y(i,2);
   rd(i) =  k2*y(i,1)-k3*y(i,2)*y(i,7);  
   re(i) =  k3*y(i,2)*y(i,7);
   Scd(i)=y(i,6)/(y(i,9)+0.00001);
   Sde(i)=y(i,9)/(y(i,10)+0.00001);
end
nexttile
plot(t,X(1:z(1,1)),t,y(:,11),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('Conversion')
leg=legend('X','Xbar','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Conversion profile')
axis([0 6 0 1])
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 6 -3 3])
grid on
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 6 0 20])
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;
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;
X=(cao-ca)/cao;
E1 =  0.47219*t^4-1.30733*t^3+0.31723*t^2+0.85688*t+0.20909;
E2 =  3.83999*t^6-58.16185*t^5+366.20970*t^4-1224.66963*t^3+2289.84857*t^2-2265.62125*t+925.46463;
E3 =  0.00410*t^4-0.07593*t^3+0.52276*t^2-1.59457*t+1.84445;
if t<=1.82
    E=E1;
elseif t<=2.8
    E=E2;
else
  E=E3;
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