% LEP 8-7 : Complex Reactions in a Semi-batch Reactor
%

tspan = [0 100]; % Range for the independent variable time
y0 = [200; 0; 0; 0]; % Initial values for the dependent variables i.e. Na,Nb,Nc,and Nd 
k1a = 10; 
k2c = 15; 
Fao = 3; 
vo=10;
Vo=1000;
[t,y]=ode45(@(t,y)ODEfun(t,y,k1a,k2c,Vo,vo,Fao),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,800,400])
nexttile
plot(t,y(:,2),t,y(:,1),t,y(:,4),t,y(:,3),'LineWidth',2)
leg=legend('N_A','N_B','N_C','N_D');
leg.ItemTokenSize=[10,10];
xlabel('t(min)')
ylabel('N_i (moles)')
axis([0 100 0 300])
title('Number of moles as a function of time')
grid on
z=size(y);
s=size(t);
Nbo=200;
% Calculating values of selectivity
for i=1:z(1,1)
    if t(i,1)>0.0001
    Scd(i)=y(i,4)/y(i,3);
    else 
        Scd(i)=0;
    end
    X(i)=1-y(i,1)/Nbo;
end
nexttile
plot(t,Scd(1:s(1,1)),'LineWidth',2)
axis([0 100 0 10^9])
xlabel('t(min)')
ylabel('Selectivity')
leg=legend('S_C_/_D');
leg.ItemTokenSize=[10,10];
title('Selectivity as a function of time')
grid on
nexttile
plot(t,X(1:s(1,1)),'LineWidth',2)
axis([0 100 0 1])
xlabel('t(min)')
ylabel('Conversion')
leg=legend('X');
leg.ItemTokenSize=[10,10];
title('Conversion as a function of time')
grid on

%

function dYfuncvecdt = ODEfun(t,Yfuncvec,k1a,k2c,Vo,vo,Fao); 
Nb = Yfuncvec(1); 
Na = Yfuncvec(2); 
Nd = Yfuncvec(3); 
Nc = Yfuncvec(4); 
% Explicit equations
V = Vo + vo * t; 
Ca = Na / V; 
Cb = Nb / V; 
r1a = 0 - (k1a * Ca * Cb ^ 2); 
Cc = Nc / V; 
r1b = 2 * r1a; 
rb = r1b; 
r2c = 0 - (k2c * Ca ^ 2 * Cc ^ 3); 
r2a = 2 / 3 * r2c; 
r2d = -1 / 3 * r2c; 
r1c = 0 - r1a; 
rd = r2d; 
ra = r1a + r2a; 
Cd = Nd / V; 
rc = r1c + r2c; 
if (t > 0.0001) 
    Scd = Nc / Nd;
else
    Scd = 0;
end 
% Differential equations
dNbdt = rb * V; 
dNadt = ra * V + Fao; 
dNddt = rd * V; 
dNcdt = rc * V; 
dYfuncvecdt = [dNbdt; dNadt; dNddt; dNcdt]; 
end

%