% Example LEP-WEB17-1a: Reactions in a Series in a PFR
%

Vspan = [0 12.6]; % Range for the independent variable
y0 = [1 0 0]; % Initial values for the dependent variables
k1=1;
k2=5;
vo=10;
Cao=1;
[v,y]=ode45(@(v,y)ODEfun(v,y,k1,k2,vo,Cao),Vspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,500])
nexttile
plot(v,y(:,1),v,y(:,2),v,y(:,3),'LineWidth',2)
grid on
xlabel('Volume ')
ylabel('Concentration')
leg=legend('C_A','C_B','C_C');
leg.ItemTokenSize=[10,10];
title('Concentration profile')
xlim([0 12.6])
ylim([0 1])
z=size(y);
 for i=1 : z(1,1)
    ra(i) =  -k1*y(i,1);
    rc(i) =  k2*y(i,2);
    rb(i) =  k1*y(i,1)-k2*y(i,2);
    X(i)=1-y(i,1)/Cao;  
 end
 nexttile
plot(v,X(1:z(1,1)),'LineWidth',2)
grid on
xlabel('Volume')
ylabel('Conversion')
leg=legend('X');
leg.ItemTokenSize=[10,10];
title('Conversion profile')
xlim([0 12.6])
ylim([0 1])
nexttile
plot(v,-ra(1:z(1,1)),v,-rb(1:z(1,1)),v,-rc(1:z(1,1)),'LineWidth',2)
grid on
xlabel('Volume')
ylabel('Rate')
leg=legend('-r_A','-r_B','-r_C');
leg.ItemTokenSize=[10,10];
title('Reaction rate profile')
xlim([0 12.6])
ylim([-2 2])

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,k1,k2,vo,Cao); 
ca= Yfuncvec(1);
cb= Yfuncvec(2);
cc= Yfuncvec(3);
% Explicit equations
ra =  -k1*ca;
rc =  k2*cb;
rb =  k1*ca-k2*cb;
% Differential equations
dcadV = (ra/vo);
dcbdV = (rb/vo);
dccdV = (rc/vo);
dYfuncvecdV=[dcadV;dcbdV;dccdV];
end