% Example 17-2: Mean Conversion Calculation in a Real Reactor
%

tspan = [0 14]; % Range for the independent variable
y0 = [ 0 0 0 0 0]; % Initial values for the dependent variables
k=0.1;
CA0=1;
n=1;
[t,y]=ode45(@(t,y)ODEfun(t,y,k,CA0,n),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,500])
nexttile
plot(t,y(:,1),t,y(:,2),t,y(:,5),'LineWidth',1.5);
xlabel('time (mins)')
ylabel('Conversion')
leg=legend('X','X_{seg}','X_{seg} for LFR RTD','Location','NorthWest');
axis([0 14 0 1])
leg.ItemTokenSize=[10,10];
title('Conversion profile')
grid on
nexttile
plot(t,y(:,3),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('t_m (mins)')
leg=legend('t_m','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('Mean residence time profile')
axis([0 14 0 6])
grid on
nexttile
plot(t,y(:,4),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('F')
leg=legend('F','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('F curve')
axis([0 14 0 1])
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,k,CA0,n); 
X= Yfuncvec(1);
Xseg= Yfuncvec(2);
tm= Yfuncvec(3);
F= Yfuncvec(4);
XsegLFR=Yfuncvec(5);
% Explicit equations
C1= 0.0038746 + 0.2739782*t + 1.574621*t^2 - 0.2550041*t^3;
C2= -33.43818 + 37.18972*t - 11.58838*t^2 + 1.695303*t^3 - 0.1298667*t^4 + 0.005028*t^5 - 7.743*10^-5*t^6;
if t <= 4 && t>=0
    C= C1;
elseif t > 4 && t<=14
     C= C2;
end
E=C/51;
if(t<5.1/2) 
  E_LFR =0;
else
   E_LFR= ((5.1^2)/(2*t^3+0.0001));
end
% Differential equations
dXdt=k*(CA0^n)*((1-X)^n)/CA0;
dXsegdt = X*E; 
dtmdt = t*E; 
dFdt = E;
dXsegLFR= X*E_LFR; 
dYfuncvecdt=[dXdt;dXsegdt;dtmdt;dFdt;dXsegLFR];
end