% Example 17-4:  Conversion Bounds for a Nonideal Reactor
%

tspan = [0 200]; % Range for the independent variable
y0 = [ 0 0 ]; % Initial values for the dependent variables
k=0.01;
Cao=8;
n=2;
[t,y]=ode45(@(t,y)ODEfun(t,y,k,Cao,n),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,800,700])
nexttile
plot(t,y(:,1),t,y(:,2),'LineWidth',1.5);
xlabel('time (mins)')
ylabel('Conversion')
leg=legend('X','X_{seg}','Location','NorthEast');
axis([0 200 0 1])
leg.ItemTokenSize=[10,10];
title('Conversion profile')
grid on
nexttile
z=size(y);
for i=1:z(1,1)
  Ca(i) = Cao*(1-y(i,1));
  ra(i) = -k*Ca(i).^2;
end
plot(t,Ca(1:z(1,1)),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('C_A (mol/dm^3)')
leg=legend('C_A','Location','NorthWest');
leg.ItemTokenSize=[10,10];
title('Concentration curve')
axis([0 200 0 8])
grid on
nexttile
plot(t,-ra(1:z(1,1)),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('-r_A (mol/dm^3.min)')
leg=legend('-r_A','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Rate curve')
axis([0 200 0 0.7])
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,k,Cao,n); 
X= Yfuncvec(1);
Xseg= Yfuncvec(2);
% Explicit equations
Ca = Cao*(1-X);
ra = -k*Ca^2;
E1 = 4.44658e-10*t^4-1.1802e-7*t^3+1.35358e-5*t^2 -0.000865652*t+0.028004;
E2 = -2.64e-9*t^3+1.3618e-6*t^2-0.00024069*t+0.015011;
if t<=70
    E=E1;
else 
    E=E2;
end
% Differential equations
dXdt=-ra/Cao;
dXsegdt = X*E; 
dYfuncvecdt=[dXdt;dXsegdt];
end