%
% LEP 1-3 b: Reactor Volume for Reversible Reaction
%

Vspan = [0 100];
y0 = [10 0]; 
k = 0.23;
v0 = 10;
Ke=2;
[v,y]=ode45(@(v,y)ODEfun(v,y,k,v0,Ke),Vspan,y0);
tiledlayout(1,2)
set (gcf,'Position',[0,0.1,800,300])
nexttile
plot(v,y(:,1),v,y(:,2),'LineWidth',2)
grid on
xlabel('Volume (dm^3)')
ylabel('Concentration (mol/dm^3)')
leg=legend('C_A','C_B');
leg.ItemTokenSize=[10,10];
title('Concentration profile')
xlim([0 100])
ylim([0 10])
z=size(y);
for i=1 : z(1,1)
    % Calculating the values of rate from its expression 
       ra(i) = -k*(y(i,1)-y(i,2)/Ke);
       rb(i) = -ra(i);
 end
nexttile
plot(v,-ra(1:z(1,1)),v,-rb(1:z(1,1)),'LineWidth',2)
grid on
xlabel('Volume (dm^3)')
ylabel('Rate (mol/dm^3/min)')
leg=legend('-r_A','-r_B');
leg.ItemTokenSize=[10,10];
title('Reaction rate profile')
xlim([0 100])
ylim([-10 10])

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,k,v0,Ke)
Ca = Yfuncvec(1); 
Cb = Yfuncvec(2);
% Explicit equations
ra = -k * (Ca-Cb/Ke);
rb = -ra;
% Differential equations
dCadV = ra / v0;
dCbdV = rb / v0;
dYfuncvecdV =[dCadV; dCbdV]; 
end