% step4s.m
% solve the nonlinear reactor equations
% simplified version

% 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)

figure(4), clf                % clear figure
plot(t,X);                    % the solution curves 

figure(5), clf 
plot(X(:,1),X(:,2));          % the phase portrait, X1 versus X2

