% LEP 8-3 : Series Reaction in a Batch Reactor
%

tspan = [0 20]; % Range for the time of the reaction 
y0 = [2; 0; 0]; % Initial values for the dependent variables i.e. Ca,Cb,Cc
k1o=0.5;
k2o=0.2;
E1=0.8;
E2=0.6;
T=450;
[w,y]=ode45(@(w,y)ODEfun(w,y,k1o,k2o,E1,E2,T),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,800,400])
nexttile
plot(w,y(:,1),w,y(:,2),w,y(:,3),'LineWidth',2)
leg=legend('C_A','C_B','C_C');
xlabel('time (hr)')
ylabel('C_i (mol/dm^3)')
title('Concentration Profiles')
leg.ItemTokenSize=[10,10];
axis([0 20 0 2])
grid on
z=size(y);
Cao=2;
% Calculating values of selectivity
for i=1:z(1,1)
     Yield(i)=y(i,2)/(Cao-y(i,1));
     X(i)=(Cao-y(i,1))/Cao;
    if w(i,1)>0.0001
    Scd(i)=y(i,2)/y(i,3);
    else 
        Scd(i)=0;
    end
end
nexttile
plot(w,Scd(1:z(1,1)),'LineWidth',2)
xlabel('time (hr)')
ylabel('Selectivity')
leg=legend('S_B_/_C');
leg.ItemTokenSize=[10,10];
axis([0 20 0 5])
title('Selectivity Profile')
grid on
nexttile
plot(w,Yield(1:z(1,1)), 'LineWidth',2)
xlabel('time (hr)')
ylabel('Yield')
leg=legend('Y_B');
axis([0 20 0 1])
title('Yield Profile')
leg.ItemTokenSize=[10,10];
grid on
nexttile
plot(w,X(1:z(1,1)), 'LineWidth',2)
xlabel('time (hr)')
ylabel('Conversion')
leg=legend('X');
axis([0 20 0 1])
title('Conversion Profile')
leg.ItemTokenSize=[10,10];
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,k1o,k2o,E1,E2,T); 
Ca = Yfuncvec(1); 
Cb = Yfuncvec(2); 
Cc = Yfuncvec(3); 
% Explicit equations
R=1.987*10^-3;
k1=k1o*exp((E1/R)*(1/300-1/T));
k2=k2o*exp((E2/R)*(1/300-1/T));
ra=-k1*Ca;
rb=k1*Ca-k2*Cb;
rc=k2*Cb;
% Differential equations
dCadt = ra; 
dCbdt = rb; 
dCcdt = rc; 
dYfuncvecdt = [dCadt; dCbdt; dCcdt]; 
end