% LEP 8-1 :Selectivity Maximization for Trambouze Reactions (PFR)
%

tau = [0 300]; % Range for the volume of the reactor 
y0 = [0.112;0.0783;0.132;0.0786]; % Initial values for the dependent variables
k1 =0.0001;
k2 = 0.0015;
k3 = 0.008;
[v y]=ode45(@(v,y)ODEfun(v,y,k1,k2,k3),tau,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,600,400])
nexttile
plot(v,y(:,1),v,y(:,2),v,y(:,3),v,y(:,4),'LineWidth',2);
leg=legend('C_A','C_X','C_B','C_Y');
leg.ItemTokenSize=[10,10];
axis([0 300 0 0.3])
xlabel('tau(s)')
ylabel('C(mol/dm^3)')
title('Concentration Profile')
grid on
z=size(y);
Cao = .4;

% Calculating values of selectivity
for i=1:z(1,1)
     Sbxy(i)=y(i,3)/(y(i,2)+y(i,4));
     X(i)=(Cao-y(i,1))/Cao;
end
nexttile
plot(v,Sbxy(1:z(1,1)),'LineWidth',2)
xlabel('tau(s)')
ylabel('Selectivity')
leg=legend('S_B_/_{XY}');
leg.ItemTokenSize=[10,10];
axis([0 300 0.4 0.9])
title('Selectivity Profile')
grid on
nexttile
plot(v,X(1:z(1,1)), 'LineWidth',2)
xlabel('tau(s)')
ylabel('Conversion')
leg=legend('X');
axis([0 300 0 1])
title('Conversion Profile')
leg.ItemTokenSize=[10,10];
grid on

%

function dYfuncvecdtau = ODEfun(tau,Yfuncvec,k1,k2,k3); 
Ca =Yfuncvec(1);
Cx =Yfuncvec(2);
Cb = Yfuncvec(3);
Cy = Yfuncvec(4);

% Differential equations

dCadtau = -k1-k2*Ca-k3*Ca^2;
dCxdtau = k1;
dCbdtau = k2*Ca;
dCydtau = k3*Ca^2;
dYfuncvecdtau=[dCadtau;dCxdtau;dCbdtau;dCydtau];
end