% Example 17-3: Mean Conversion in a Laminar flow Reactor
%

tspan = [0 20000]; % Range for the independent variable
y0 = [ 0 0 ]; % Initial values for the dependent variables
k=0.00493;
Cao=0.75;
n=2;
tau=1000;
[t,y]=ode45(@(t,y)ODEfun(t,y,k,Cao,n,tau),tspan,y0);
tiledlayout(1,2)
set (gcf,'Position',[0,0.1,1000,300])
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 20000 0 1])
leg.ItemTokenSize=[10,10];
title('Conversion profile')
grid on
nexttile
z=size(y);
t1=tau/2;
for i=1:z(1,1)
   E2(i) = tau^2/2/(t(i).^3+.00001);
if(t(i)<t1) 
  E(i) =0;
else
   E(i)= E2(i);
end 
end
plot(t,E(1:z(1,1)),'LineWidth',1.5)
xlabel('time (mins)')
ylabel('E')
leg=legend('E','Location','NorthEast');
leg.ItemTokenSize=[10,10];
title('E curve')
axis([0 20000 0 0.005])
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,k,CA0,n,tau); 
X= Yfuncvec(1);
Xbar= Yfuncvec(2);
% Explicit equations
t1=tau/2;
E2 = tau^2/2/(t^3+.00001);
if(t<t1) 
  E =0;
else
   E= E2;
end
% Differential equations
dXdt=k*(CA0^n)*((1-X)^n)/CA0;
dXbardt = X*E; 
dYfuncvecdt=[dXdt;dXbardt];
end