% Example LEP-17-6b: Maximum Mixedness model with Bimodal distribution (Multiple Reactions)
%

zspan = [0 6]; % Range for the independent variable
y0 = [1 1 0 0.99 0 0 ]; % 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(:,5),z,y(:,6),'LineWidth',1.5);
xlabel('z(mins)')
ylabel('Concentration')
leg=legend('C_A ','C_B','C_C','C_D','C_E','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Concentration profile')
axis([0 6 0 1])
grid on
a=size(y);
cao=1;
for i=1:a(1,1)
    X(i)=(cao-y(i,1))/cao;
    Scd(i)=y(i,3)/y(i,5);
    Sde(i)=y(i,3)/y(i,6);
    rb(i) = -k1*y(i,1)*y(i,2)-k3*y(i,2)*y(i,5);
    ra(i) = -k1*y(i,1)*y(i,2)-k2*y(i,1);
    rd(i) = k2*y(i,1)-k3*y(i,2)*y(i,5);
    rc(i) = k1*y(i,1)*y(i,2);
    re(i) = k3*y(i,2)*y(i,5);
end
nexttile
plot(z,X(1:a(1,1)),'LineWidth',1.5)
xlabel('z (mins)')
ylabel('$\bar{X}$','Interpreter','Latex')
leg=legend('Xbar ','Location','NorthWest');
leg.ItemTokenSize=[10,10];
title('Conversion profile')
axis([0 6 0 1])
grid on
nexttile
plot(z,Scd(1:a(1,1)),z,Sde(1:a(1,1)),'LineWidth',1.5)
xlabel('z (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
nexttile
plot(z,-ra(1:a(1,1)),z,-rb(1:a(1,1)),z,-rc(1:a(1,1)),z,-rd(1:a(1,1)),z,-re(1:a(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 6 -3 3])
grid on

%

function dYfuncvecdz = ODEfun(z,Yfuncvec,k1,k2,k3); 
ca= Yfuncvec(1);
cb= Yfuncvec(2);
cc= Yfuncvec(3);
F= Yfuncvec(4);
cd= Yfuncvec(5);
ce= Yfuncvec(6);
% Explicit equations
cbo = 1;
cao = 1; 
cco = 0; 
cdo = 0; 
ceo = 0; 
lam = 6-z;
rc = k1*ca*cb;
re = k3*cb*cd;
E1 = 0.47219*lam^4-1.30733*lam^3+0.31723*lam^2+0.85688*lam+0.20909;
E2 = 3.83999*lam^6-58.16185*lam^5+366.2097*lam^4-1224.66963*lam^3+2289.84857*lam^2-2265.62125*lam+925.46463;
E3 = 0.00410*lam^4-0.07593*lam^3+0.52276*lam^2-1.59457*lam+1.84445;
rb = -k1*ca*cb-k3*cb*cd;
ra = -k1*ca*cb-k2*ca;
rd = k2*ca-k3*cb*cd;
if lam<=1.82
E = E1;
elseif lam<=2.8
    E=E2;
else
 E=E3;
end
EF = E/(1-F);
% Differential equations
dcadz = -(-ra+(ca-cao)*EF);
dcbdz = -(-rb+(cb-cbo)*EF);
dccdz = -(-rc+(cc-cco)*EF); 
dFdz = -E;
dcddz = -(-rd+(cd-cdo)*EF);
dcedz = -(-re+(ce-ceo)*EF);
dYfuncvecdz=[dcadz;dcbdz;dccdz;dFdz;dcddz;dcedz];
end