% step4.m
% solve the nonlinear reactor equations

% first check the eigenvalues
beta
arr=arrhen(Xbar(2)); 
arrprime=arr*gamma*Xbar(2)^(-2);
A=[-1-arr -Xbar(1)*arrprime; alpha*arr -1-beta+alpha*Xbar(1)*arrprime];
lambda=eig(A)

% now solve the equations
Ubar=Xbar(2)-(1-Xbar(2)+alpha*(1-Xbar(1)))/beta 
x0, t1, u0 
X0=Xbar+x0
[t,X]=ode45('reactor',[0 t1],X0);
% [t,X]=ode15s('reactor',[0 t1],X0);  % if ode45 is slow (stiff system)

% plot the solutions 
figure(4), clf  % clear figure
subplot(2,1,1)
% the stationary solution, two red lines
plot(t,Xbar(1)*ones(size(t)),'r',t,Xbar(2)*ones(size(t)),'r');
hold on; 
plot(t,X);      % the solution curves as functions of time

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

subplot(2,1,2) 
plot(Xbar(1),Xbar(2),'xr');   % the stationary solution, a red cross
hold on;
plot(X(:,1),X(:,2));          % the phase plane, X1 versus X2
% axis ([0 1 1 1.2]) 
% axis square
title([ 'x marks the stationary solution' ])

% the following commands plot the solutions of the nonlinear 
% and linearized systems in the same figure 
% figure(5), clf 
% plot(X(:,1),X(:,2),'b',Xbar(1),Xbar(2),'xr', ... 
%    Xbar(1)+x(:,1),Xbar(2)+x(:,2),'g')  
% axis ([0 1 1 1.2]) 
