% LEP 13-1a : Adiabatic Batch Reactor (Part a)
%

tspan = [0 2500]; % Range for the independent variable
y0 = [0]; % Initial values for the dependent variables X
k1a=0.000273;
Ea=18000.23;
To=286;
dH=-20202;
NA0=54.8;
NB0=555;
NM=98.8;
R=1.987;
[t y]=ode45(@(t,y)ODEfun(t,y,k1a,Ea,To,dH,NA0,NB0,NM,R),tspan,y0);
tiledlayout(1,2)
set (gcf,'Position',[0,0.1,900,300])
nexttile
plot(t,y(:,1),'LineWidth',2);
leg=legend('X');
ylabel('Conversion');
xlabel('t(s)');
axis([0 2500 0 1])
title('Adiabatic Operation: Batch Reactor')
leg.ItemTokenSize=[10,10];
grid on
z=size(y);
thetaB=NB0/NA0;
thetaM=NM/NA0;
Cps=35+18*thetaB+19.5*thetaM;
% Calculating values of selectivity
for i=1:z(1,1)
     T(i)=To+(-dH/Cps)*y(i);
end
nexttile
plot(t,T(1:z(1,1)),'LineWidth',2)
leg=legend('T');
leg.ItemTokenSize=[10,10];
title('Adiabatic Operation: Batch Reactor')
ylabel('T(K)');
xlabel('t(s)');
axis([0 2500 280 480]);
grid on

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,k1a,Ea,To,dH,NA0,NB0,NM,R); 

thetaB=NB0/NA0;
thetaM=NM/NA0;
Cps=35+18*thetaB+19.5*thetaM;
X = Yfuncvec(1);
% Explicit equations
T=To+(-dH/Cps)*X;
k = k1a * exp((Ea/R) * (T - 297) / (T * 297)); 
% Differential equations
dXdt = k*(1-X); 
dYfuncvecdV = [dXdt];  
end