% LEP 10-4 : Calculating Conversion with Catalyst Decay in Batch Reactors
%

tspan = [0 500]; % Range for the independent variable t
y0 = [0,0]; % Initial values for the dependent variables i.e X and Xd
k1=0.01;
k1d=0.1;
Ea=20000;
Ed=75000;
T=293;
[t,y]=ode45(@(t,y)ODEfun(t,y,k1,k1d,Ea,Ed,T),tspan,y0);
plot(t,y(:,1),t,y(:,2),'LineWidth',2)
ylabel('Conversion');
xlabel('time (s)');
leg=legend('X','X_d');
leg.ItemTokenSize=[10,10];
axis([0 500 0 1])
title('Conversion with Catalyst Decay in Batch Reactors')
grid on

%

function dYfuncvecdt = ODEfun(t,Y,k1,k1d,Ea,Ed,T); 
X = Y(1);
Xd=Y(2);
R=1.987;
T1=300;
T2=300;
% Explicit equations
k = k1*exp((Ea/R)*((1/T1) - (1/T)));
kd = k1d*exp((Ed/R)*((1/T2) - (1/T)));
% Differential equations
dXdt= k*(1-X); 
dXddt=k*(1-Xd)/(1+kd*t);
dYfuncvecdt = [dXdt;dXddt]; 
end
