% Synthron Explosion
%

tspan = [0 3000]; % Range for the independent variable i.e.time
Cao = 2.22 ;
Cco = 3.52 ;
Cto = 2.92;
Vo = 0.73;
vo=0.00015;
Cai=2.12;
Nci=2.64;
Nti=2.2;
UA=5;
Cpa = 245.75;
Cpc = 161.3;
Cpt = 165.60;
Cpan = 231;
y0 = [Cai; Nci;Nti;0;355]; % Initial values for the dependent variables 
tiledlayout(2,2)
set (gcf,'Position',[0,0.1,800,600])
nexttile
[t y]=ode45(@(t,y)ODEfun(t,y,Cao,Cco,Cto,Vo,vo,UA,Cpa,Cpc,Cpt,Cpan),tspan,y0);
plot(t,y(:,5),'LineWidth',1.5);
ylabel('T(K)');
xlabel('t(s)');
axis([0 3000 300 600]);
title('Temperature profile')
grid on
nexttile
plot(t,y(:,1),t,y(:,4),'LineWidth',2);
leg=legend('C_A','C_{An}');
leg.ItemTokenSize=[10,10];
title('Concentration profile')
ylabel('Concentration (kmol/m^3)');
xlabel('t(s)');
axis([0 3000 0 3]);
grid on
z=size(t);
for i=1:z(1,1)
    V(i) = Vo+(vo*t(i));
    k(i) = 4.01*10^3*exp(-29000/(8.314*y(i,5)));
    ra(i) = -k(i).*(y(i,1).^2);
    Qgs(i) = -ra(i).*V(i)*64512; 
end
nexttile
plot(t,Qgs(1:z(1,1)),'LineWidth',2);
leg=legend('Q_g');
leg.ItemTokenSize=[10,10];
title('Heat-generated profile')
ylabel('Q(J/s)');
xlabel('t(s)');
axis([0 3000 0 400000]);
grid on

%

function dYfuncvecdV = ODEfun(t,Yfuncvec,Cao,Cco,Cto,Vo,vo,UA,Cpa,Cpc,Cpt,Cpan); 
Ca = Yfuncvec(1); 
Nc = Yfuncvec(2); 
Nt = Yfuncvec(3);
Can = Yfuncvec(4);
T = Yfuncvec(5);
% Explicit equations
V = Vo+(vo*t);
Na= Ca*V;
Nan = Can*V;
Fao = Cao*vo;
Fco = Cco*vo;
Fto= Cto*vo;
k = 4.01*10^3*exp(-29000/(8.314*T));
ra = -k*(Ca^2);
Qgs = -ra*V*64512;
Qrs1 = ((Fao*Cpa)+(Fco*Cpc)+(Fto*Cpt))*(T-310);
Qrs2 = UA*(T-298);
Qrs= Qrs1+Qrs2+8.3;
NCp = (Na*Cpa)+(Nc*Cpc)+(Nt*Cpt)+(Nan*Cpan);

% Differential equations
dCadt = (vo/V)*(Cao-Ca) + ra;
dNcdt = Cco*vo;
dNtdt = Cto*vo;
dCandt = -ra - (Can*(vo/V));
dTdt = (Qgs-Qrs)/NCp;
dYfuncvecdV = [dCadt; dNcdt; dNtdt;dCandt;dTdt];  
end