function SemiBatchDispl
clear, clc, format short g, format compact
tspan = [0.0001 7.2E+04]; % Range for the independent variable 
y0 = [0; 0; 260.; 273.15]; % Initial values for the dependent variables 
disp(' Variable values at the initial point ');
disp([' t    = ' num2str(tspan(1))]);
disp('           y                  dy/dt         ');
disp([y0 ODEfun(tspan(1),y0)]);
[t,y]=ode45(@ODEfun,tspan,y0);
for i=1:size(y,2)
    disp([' Solution for dependent variable y' int2str(i)]);
    disp(['              t                  y' int2str(i)]);
    disp([t y(:,i)]);
    plot(t,y(:,i));
    title([' Plot of dependent variable y' int2str(i)]);
    xlabel(' Independent variable (t)');
    ylabel([' Dependent variable y' int2str(i)]);
    pause
end
%- - - - - - - - - - - - - - - - - - - - - -
function dYfuncvecdt = ODEfun(t,Yfuncvec); 
Np = Yfuncvec(1); 
Nx = Yfuncvec(2); 
Tr = Yfuncvec(3); 
Tcool = Yfuncvec(4); 
%Initial concentr. of nitrosonium ion Y=Nb0/NaF 
Y = 0.035; 
%Molar mass of HNO3 [kg/kmol] 
Mw = 63; 
%Initial volume in a reactor [m3] 
Vr0 = 1.5; 
%Final volume of the dose [m3] 
Vdos1 = 0.6; 
%Density of pure nitric acid [kg/m3] 
RhoAcid = 1500; 
%Initial mass concentr of nitr. acid sol. [%] 
Percent = 0.6; 
%Molar mass of 2-octanol [kg/kmol] 
MwOctan = 130.23; 
%dosing time [s], 10h 
tdos = 36000; 
%Density of 2-octanol [kg/m3] 
RhoOctan = 820.7; 
%Volume fraction of dispersed phase 
Epsd = Vdos1 / (Vdos1 + Vr0); 
%Pre-exponential factor reaction 1 [m3/kmol/s] 
maA1 = 10 ^ 5; 
%Pre-exponential factor reaction 2[m3/kmol/s] 
mpA2 = 10 ^ 10; 
%Activation temperature reaction 1 [K] 
E1perR = 11300; 
%Activation tempetature reaction 2 [K] 
E2perR = 12000; 
%Hammett's reaction rate coeff. reaction 1 
m1 = 6.6; 
%Hammett's reaction rate coeff. reaction 2 
m2 = 2.2; 
%Initial number of mole of HNO3 [kmole] 
NnO = Vr0 * Percent * RhoAcid / Mw; 
%Dimensionless time up to t=tdos 
if (t <= tdos) 
    Theta = (t / tdos) ; 
else
    Theta = (1); 
end 
%Total amount of 2-octanol (a) fed [kmol] 
NaF = Vdos1 * RhoOctan / MwOctan; 
%Concentr. of HNO3 in the aq. phase [kmol/m3] 
CnAq = (NnO - Y * NaF - Np - 2 * Nx) / Vr0; 
%Number of moles of HNO3 [kmol] 
Nn = CnAq * Vr0; 
%Concentr. of (P) in org phase [kmol/m3] 
CpOrg = Np / (Vdos1 * Theta); 
%Concentr of a in org phase [kmole/m3] 
CaOrg = (Theta * NaF - Np - Nx) / (Vdos1 * Theta); 
%Concentr. of (B) in aq. phase [kmole/m3] 
CbAq = (Np + Y * NaF) / Vr0; 
%Mass concentr. of nitric acid sol [%/100%] 
wt = Nn * Mw / (Vr0 * RhoAcid); 
%Hammett's acidity function 
H = -.6221 - 3.7214 * wt - 1.5714 * wt ^ 2; 
%Specific reaction rate 1 
k1 = maA1 * exp(-E1perR / Tr - m1 * H); 
%Specific reaction rate 2 
k2 = mpA2 * exp(-E2perR / Tr - m2 * H); 
%Reaction rate of a and b to p[kmol/m3/s] 
r1 = k1 * CaOrg * CbAq * (1 - Epsd); 
%Reaction rate of p and b to x[kmol/m3/s] 
r2 = k2 * CpOrg * CbAq * (1 - Epsd); 
%Specific heat of reaction 2 [J/kmole] 
Hnone = 520 * 10 ^ 6; 
%Temperature of feed dose [K] 
Tdos = 293.15; 
%Initial cool. surface heat trans. coeff.[W/K] 
UA0 = 1500; 
%Initial heat capacity of the system [J/K] 
Gamma0 = 5.4 * 10 ^ 6; 
%Specific heat of reaction 1 [J/kmol] 
Hnol = 160 * 10 ^ 6; 
%Heat of reaction, 2 [W] 
Qnone = r2 * Vr0 * Hnone / (1 - Epsd); 
%Volumetric flow rate of the feed [m3/s] 
Phi = Vdos1 / tdos; 
%Heat capacity of dose [J/m3/K] 
RhoCPdos = 2 * 10 ^ 6; 
%Heat input due to reactant addition [W] 
Qdos = Phi * RhoCPdos * (Tdos - Tr); 
%Final cool. surface heat trans. coeff. [W/K] 
UA1 = 2100; 
%Total heat capacity of the system [J/K] 
Gamma = Gamma0 + RhoCPdos * Phi * t; 
%Heat of reaction, 1 [W] 
Qnol = r1 * Vr0 * Hnol / (1 - Epsd); 
%Sum of the heat of reaction the reactions [W) 
Qr = Qnol + Qnone; 
%Cooling surface heat transfer coefficient [W/K] 
UAcool = UA0 + (UA1 - UA0) * Theta; 
%Heat removed by the cooling jacket [W] 
Qcool = UAcool * (Tcool - Tr); 
%Flow rate of coolant [m3/s] 
Fw = 100 / 60 * 10 ^ (-3); 
%Initial coolant temperature [K] 
Tcool_IN = 260; 
%The density of coolant [kg/m3] 
RhoCoolant = 1000; 
%Heat capacity of coolant [J/kg/K] 
CpCoolant = 4183; 
%Volume of the jacket [m3] 
Vj = 1.5; 
%Mole balance for 2-octanone (P) 
dNpdt = (r1 - r2) * Vr0 / (1 - Epsd); 
%Mole balance for carboxylic acids (X) 
dNxdt = r2 * Vr0 / (1 - Epsd); 
%Reactor energy balance (Tr in K) 
dTrdt = (Qr + Qdos + Qcool) / Gamma; 
%Jacket energy balance (T in K) 
dTcooldt = (Fw * (Tcool_IN - Tcool) - Qcool / (RhoCoolant * CpCoolant)) / Vj; 
dYfuncvecdt = [dNpdt; dNxdt; dTrdt; dTcooldt];