% Example 17-5:\,\textrm{Using}\,\textrm{Software}\,\textrm{to}\,\textrm{Make}\,\textrm{Maximum}\,\textrm{Mixedness}\,\textrm{Model}\,\textrm{Calculations}
%

zspan = [0 200]; % Range for the independent variable
y0 = [0] ; % Initial values for the dependent variables
k=0.01;
Cao=8;
n=2;
[z,y]=ode45(@(z,y)ODEfun(z,y,k,Cao,n),zspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,500])
nexttile
plot(z,y(:,1),'LineWidth',1.5);
xlabel('z (mins)')
ylabel('Conversion')
leg=legend('X');
axis([0 200 0 1])
grid on
leg.ItemTokenSize=[10,10];
title('conversion profile')
p=size(y);
for i=1:p(1,1)
   Ca(i)=Cao*(1-y(i,1));  
   ra(i) =  -k*Ca(i)^n;

end
nexttile
plot(z,Ca(1:p(1,1)),'LineWidth',1.5)
xlabel('z(mins)')
ylabel('Concentration')
leg=legend('C_A (mol/dm^3)','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Concentration profile')
axis([0 200 0 8])
grid on
nexttile
plot(z,-ra(1:p(1,1)),'LineWidth',1.5)
xlabel('z(mins)')
ylabel('Rate (mol/dm^3.min)')
leg=legend('-r_A','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Rate profile')
axis([0 200 0 0.8])
grid on

%

function dYfuncvecdz = ODEfun(z,Yfuncvec,k,Cao,n); 
X= Yfuncvec(1);
% Explicit equations
Ca=Cao*(1-X);
ra= -k*Ca^n;
lam = 200 - z;
E1 = 4.44658e-10*lam^4-1.1802e-7*lam^3+1.35358e-5*lam^2-.000865652*lam+.028004;
E2 = -2.64e-9*lam^3+1.3618e-6*lam^2-.00024069*lam+.015011;
F1 = 4.44658e-10/5*lam^5-1.1802e-7/4*lam^4+1.35358e-5/3*lam^3-.000865652/2*lam^2+.028004*lam;
F2 = -(-9.30769e-8*lam^3+5.02846e-5*lam^2-.00941*lam+.618231-1);
if lam<=70
    E=E1;
    F=F1;
else
  E=E2;
  F=F2;
end 
EF=E/(1-F);
% Differential equations
dXdz = -(ra/Cao+EF*X);
dYfuncvecdz=[dXdz];
end