% Table T12-2: Computer Experiment
%

Vspan = [0 4500]; % Range for the independent variable
y0 = [320.; 1; 330; 0]; % Initial values for the dependent variables
alpha=0.0002;
Uarho=0.5;
mc=2000;
Cpmc=18;
Hr=-20000;
Fao=5;
thetaI=1;
CpI=40;
CpA=20;
thetaB=1;
CpB=20;
Cto=0.3;
Ea=25000;
[v y]=ode45(@(v,y)ODEfun(v,y,alpha,Uarho,mc,Cpmc,Hr,Fao,thetaI,CpI,CpA,thetaB,CpB,Cto,Ea),Vspan,y0);
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,600,400])
nexttile
plot(v,y(:,3),v,y(:,1),'LineWidth',1.5);
leg=legend('T','T_a');
leg.ItemTokenSize=[10,10];
ylabel('T(K)');
xlabel('W(kg)');
axis([0 4500 300 420]);
title('Temperature profile')
grid on
yao = 1/(1 + thetaB + thetaI);
To = 330;
Cao = yao*Cto;
sumcp = (thetaI*CpI + CpA + thetaB*CpB);
z=size(y);
for i=1:z(1,1)
Kc(i)= 1000*(exp(Hr/1.987*(1/303 - 1/y(i,3))));
ka(i)= .004*exp(Ea/1.987*(1/310 - 1/y(i,3)));
Xe(i)= ((thetaB + 1)*Kc(i)- (((thetaB + 1)*Kc(i)).^2 -4*(Kc(i) - 4)*(Kc(i)*thetaB)).^0.5)/(2*(Kc(i)- 4));
Ca(i)= Cao*(1 - y(i,4)).*y(i,2)*To./y(i,3);
Cb(i)= Cao*(thetaB -y(i,4)).*y(i,2)*To./y(i,3);
Cc(i)= Cao*2*y(i,4).*y(i,2)*To./y(i,3);
ra(i)= -ka(i)*(Ca(i)*Cb(i) - (Cc(i))^2/Kc(i));
Qg(i) = (-ra(i))*(-Hr); 
Qr(i) = Uarho*(y(i,3) - y(i,1)); 
end
nexttile
plot(v,y(:,4),v,Xe(1:z(1,1)),'LineWidth',1.5);
leg=legend('X','X_e');
leg.ItemTokenSize=[10,10];
ylabel('Conversion');
xlabel('W(kg)');
axis([0 4500 0 1]);
title('Conversion profile')
grid on
nexttile
plot(v,Qg(1:z(1,1)),v,Qr(1:z(1,1)),'LineWidth',1.5);
leg=legend('Q_g','Q_r');
leg.ItemTokenSize=[10,10];
ylabel('Q (cal/kg.s)');
xlabel('W(kg)');
axis([0 4500 0 200]);
title('Heat profile')
grid on
nexttile
plot(v,-ra(1:z(1,1)),'LineWidth',1.5);
leg=legend('-r_A');
leg.ItemTokenSize=[10,10];
ylabel('Rate (mol/kg.s)');
xlabel('W(kg)');
axis([0 4500 0 0.02]);
title('Rate profile')
grid on

%

function dYfuncvecdV = ODEfun(V,Yfuncvec,alpha,Uarho,mc,Cpmc,Hr,Fao,thetaI,CpI,CpA,thetaB,CpB,Cto,Ea); 
Ta = Yfuncvec(1); 
p = Yfuncvec(2); 
T = Yfuncvec(3); 
X = Yfuncvec(4); 
% Explicit equations
Kc= 1000*(exp(Hr/1.987*(1/303 - 1/T)));
ka= .004*exp(Ea/1.987*(1/310 - 1/T));
yao = 1/(1 + thetaB + thetaI);
To = 330;
Xe= ((thetaB + 1)*Kc- (((thetaB + 1)*Kc)^2 -4*(Kc - 4)*(Kc*thetaB))^0.5)/(2*(Kc- 4));
Cao = yao*Cto;
sumcp = (thetaI*CpI + CpA + thetaB*CpB);
Ca= Cao*(1 - X)*p*To/T;
Cb= Cao*(thetaB - X)*p*To/T;
Cc= Cao*2*X*p*To/T;
ra= -ka*(Ca*Cb - (Cc)^2/Kc);
Qg = (-ra)*(-Hr); 
Qr = Uarho*(T - Ta); 
% Differential equations
dTadW = Uarho*(T-Ta)/(mc*Cpmc); 
dpdW = -alpha/2*(T/To)/p; 
dTdW = (Qg-Qr)/(Fao*sumcp); 
dXdW = -ra/Fao; 
dYfuncvecdV = [dTadW; dpdW; dTdW; dXdW]; 
end