% Example LEP-WEB17-1f: Maximum Mixedness Model with Bimodal distribution
%

zspan = [0 6]; % Range for the independent variable
y0 = [1 0 0 0.99]; % Initial values for the dependent variables
k1=1;
k2=1;
[z,y]=ode45(@(z,y)ODEfun(z,y,k1,k2),zspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,500])
nexttile
plot(z,y(:,1),z,y(:,2),z,y(:,3),'LineWidth',1.5);
xlabel('z (mins)')
ylabel('Concentration')
leg=legend('C_Abar','C_Bbar','C_Cbar','Location','NorthEast');
axis([0 6 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);
   rb(i) =  k1*y(i,1)-k2*y(i,2);
   rc(i) =  k2*y(i,2);
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 6 0 1])
grid on
nexttile
plot(z,-ra(1:p(1,1)),z,-rb(1:p(1,1)),z,-rc(1:p(1,1)),'LineWidth',1.5)
xlabel('z(mins)')
ylabel('Rate')
leg=legend('-r_A','-r_B','-r_C','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Rate profile')
axis([0 6 -3 3])
grid on

%

function dYfuncvecdz = ODEfun(z,Yfuncvec,k1,k2); 
Ca= Yfuncvec(1);
Cb= Yfuncvec(2);
Cc= Yfuncvec(3);
F= Yfuncvec(4);
% Explicit equations
ra= -k1*Ca;
rb= k1*Ca- k2*Cb;
rc= k2*Cb;
Cao = 1;
lam = 6 - z;
Cbo = 0;
Cco = 0;
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.4646;  
E3 =  0.00410*lam^4-0.07593*lam^3+0.52276*lam^2-1.59457*lam+1.84445;  
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;
dYfuncvecdz=[dCadz;dCbdz;dCcdz;dFdz];
end