The following m-files are used in exercise on reactor stability. 

-----------------------------------------------------------------------
% data.m 
%
% TMA680 Tillampad matematik K 1996
% 1996-03-08 
% Stig Larsson

% data for exercise on reactor stability 
disp('data for exercise on reactor stability')
% given data 
disp('given data')
q =         2.79 *10^(-4)         % m^3/s 
vol =       1.20                  % m^3 
cf =        5.00 *10^3            % mol/m^3 
tf =      343.00                  % K 
deltah =  -83.70 *10^3            % J/mol 
cp =        4.19 *10^3            % J/(kg K) 
rho =       1.00 *10^3            % kg/m^3
kappa =     0.58 *10^3            % J/(m^2 s K) 
tau =  vol/q                      % s 

% measured reaction rates  
disp('measured reaction rates')
temperatures = [343 353 363 373 383 393 403]'                 % K 
rates        = 10^(-5)*[2.8 5.6 11.2 22.4 44.8 89.6 179.2]'   % s^(-1) 

% derived dimensionless variables 
disp('derived dimensionless variables') 
global temp rate; 
temp=temperatures/tf
rate=rates*tau 

-----------------------------------------------------------------------
% step1.m 
% 
global alpha beta gamma delta X1 X2 Ubar u0 A B; 

% read given data
data

% fit the rate function to measured data 
xmin=fmins('fel',[1 1]);
delta=xmin(1)
gamma=xmin(2)
plot(temp,rate,'co',temp,arrhen(temp))

% compute variables 
alpha=-deltah*cf/(rho*cp*tf)
X1=0.5
delta1=delta*X1/(1-X1);
X2=1-log(delta1)/(gamma+log(delta1)) 

-------------------------------------------------------------------
% step2.m 
% compute the eigenvalues of the Jacobian 
beta=input('give beta: ');
arr=arrhen(X2); arrprime=arr*gamma*X2^(-2);
A=[-1-arr -X1*arrprime; alpha*arr -1-beta+alpha*X1*arrprime] 
[vv,lambda]=eig(A)

% compute the corresponding value of Ubar 
Ubar=X2-(1-X2+alpha*X1*arrhen(X2))/beta

--------------------------------------------------------------------
% step3.m
% solve the linearized system
beta
B=[0;beta]; 
arr=arrhen(X2); arrprime=arr*gamma*X2^(-2);
A=[-1-arr -X1*arrprime; alpha*arr -1-beta+alpha*X1*arrprime] 
[vv,lambda]=eig(A)
Ubar=X2-(1-X2+alpha*X1*arr)/beta 
x0, t1, u0 
[t,x]=ode45('linear',0,t1,x0);

subplot(2,1,1), plot(t,x);

title(['linjarized,  beta=',num2str(beta),',  u0=',num2str(u0),',  x0=[',num2str(x0(1)),'; ',num2str(x0(2)),']'])


subplot(2,1,2), plot(x(:,1),x(:,2)); 

---------------------------------------------------------------------------
% step4.m
% solve the nonlinear system
beta
arr=arrhen(X2); arrprime=arr*gamma*X2^(-2);
A=[-1-arr -X1*arrprime; alpha*arr -1-beta+alpha*X1*arrprime] 
[vv,lambda]=eig(A)
Ubar=X2-(1-X2+alpha*(1-X1))/beta 
x0, t1, u0 
X0=[X1;X2]+x0
[t,X]=ode45('reactor',0,t1,X0);


subplot(2,1,1), plot(t,X1*ones(size(t)),'c.',t,X2*ones(size(t)),'c.',t,X);

title(['nonlinear,  beta=',num2str(beta),',  u0=',num2str(u0),',  x0=[',num2str(x0(1)),'; ',num2str(x0(2)),']'])

subplot(2,1,2), plot(X(:,1),X(:,2)); 

-----------------------------------------------------------------------------
% step5.m
% compute the state of the reactor in physical variables 
cbar=X1*cf
tbar=X2*tf-273
tkbar=Ubar*tf-273
deltatemp=tbar-tkbar
areakrit=0.898*rho*cp*vol/(kappa*tau) 
area=beta*rho*cp*vol/(kappa*tau) 
coolingpower=kappa*area*(tbar-tkbar)

aarea=1.5 
bbeta=aarea*kappa*tau/(rho*cp*vol)

-----------------------------------------------------------------------
function y=arrhen(x);
% the Arrhenius rate function
global gamma delta;
y=delta*exp(gamma-gamma./x);
end

----------------------------------------------------------------------
function y=fel(x);
% residual function to be used for fitting of 
% the Arrhenius rate function to measured data
global rate temp;
a=rate-x(1)*exp(x(2)-x(2)./temp); 
y=a'*a;
end

----------------------------------------------------------------------
function xprime=linear(t,x);
global A B; 
xprime=A*x+B*u(t); 
end 

-----------------------------------------------------------------------
function xprime=reactor(t,x);
% Ubar is the equilibrium control
% u(t) is the perturbation of the control 
% Ubar+u(t) is the control 
global alpha beta Ubar;
xprime(1)=1-x(1)      -x(1)*arrhen(x(2)); 
xprime(2)=1-x(2)+alpha*x(1)*arrhen(x(2))-beta*(x(2)-Ubar-u(t)); 
end 

----------------------------------------------------------------------
function y=u(t);
global u0;
b=1;
y=u0*cos(b*t);
end 

----------------------------------------------------------------------