% LEP P8-G1 : Oscillating Reactions
%

tspan = [0 200]; % Range for the independent variable
y0 = [0.01; 0.000025; 0.00001; 0]; % Initial values for the dependent variables i.e. Fa,Fb,Fc,Fd,and y
k0 = 0.001; 
ku = 0.01; 
k1 = 2500000000; 
k2 = 1; 
[t,y]=ode45(@(t,y)ODEfun(t,y,k0,ku,k1,k2),tspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,700,400])
nexttile
plot(t,y(:,2),t,y(:,3),'LineWidth',2)
leg=legend('C_A','C_B');
leg.ItemTokenSize=[10,10];
axis([0 200 0 0.0002])
set(gca,'ytick',linspace(0,0.0002,6))
xlabel('time (mins)')
ylabel('C_i (mol/dm^3)')
title('Concentration Profiles')
grid on
nexttile
plot(t,y(:,1),t,y(:,4),'LineWidth',2)
xlabel('time (mins)')
ylabel('C_i (mol/dm^3)')
title('Concentration Profiles')
leg=legend('C_P','C_c');
leg.ItemTokenSize=[10,10];
axis([0 200 0 0.01])
set(gca,'ytick',linspace(0,0.01,6))
grid on
Cpo=0.01;
z=size(y);
for i=1:z(1,1)
   Xp(i)=(Cpo-y(i,1))/Cpo;
end
nexttile
plot(t,Xp(1:z(1,1)),'LineWidth',2)
xlabel('time (mins)')
ylabel('Conversion')
leg=legend('X_P');
leg.ItemTokenSize=[10,10];
axis([0 200 0 1])
set(gca,'ytick',linspace(0,1,6))
title('Conversion Profile')
grid on


%

function dYfuncvecdt = ODEfun(t,Yfuncvec,k0,ku,k1,k2); 
Cp = Yfuncvec(1); 
Ca = Yfuncvec(2); 
Cb = Yfuncvec(3); 
Cc = Yfuncvec(4);  
% Explicit equations
rP=-k0*Cp;
rA= k0*Cp-k1*Ca*Cb^2-ku*Ca;
rB= ku*Ca+k1*Ca*Cb^2-k2*Cb;
% Differential equations
dCpdt = rP; 
dCadt = rA; 
dCbdt = rB; 
dCcdt = k2*Cb;  
dYfuncvecdt = [dCpdt; dCadt; dCbdt; dCcdt]; 
end

%