% Copyright (C) 2001, James B. Rawlings and John G. Ekerdt%% This program is free software; you can redistribute it and/or% modify it under the terms of the GNU General Public License as% published by the Free Software Foundation; either version 2, or (at% your option) any later version.%% This program is distributed in the hope that it will be useful, but% WITHOUT ANY WARRANTY; without even the implied warranty of% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU% General Public License for more details.%% You should have received a copy of the GNU General Public License% along with this program; see the file COPYING. If not, write to% the Free Software Foundation, 59 Temple Place - Suite 330, Boston,% MA 02111-1307, USA.globalk_mT_mEc_Afc_AC_psT_fDeltaH_RUT_a%% multiplicity parameters%% units are kmol, min, kJ, K, m^3%%k_m=0.001;T_m=298;E=8000;c_Af=2;C_p=4;rho=1000;C_ps=rho*C_p;T_f=298;T_a=T_f;DeltaH_R=-3e5;U=0;x0=[1;T_f];nc_As=250;cleartmp_table;c_Avect=linspace(0.999*c_Af,.003*c_Af,nc_As);fori=1:nc_Asc_A=c_Avect(i);opts=optimset('MaxFunEvals',2000*numel(x0),...'MaxIter',500*numel(x0));[x,fval,info]=fsolve('st_st_cA',x0,opts);theta=x(1);T=x(2);conv=(c_Af-c_A)/c_Af;%% compute the eigenvalues of the Jacobian%k=k_m*exp(-E*(1/T-1/T_m));Jac=[-1/theta-k,-k*c_A*E/(T*T);-k*DeltaH_R/C_ps,-U/C_ps-1/theta-k*c_A*DeltaH_R/C_ps*E/(T*T)];lambda=eig(Jac);a=real(lambda(1));b=imag(lambda(1));c=real(lambda(2));d=imag(lambda(2));if(a>=c)lamrow=[abcd];elselamrow=[cdab];endtmp_table(i,:)=[theta,T,conv,lamrow,info];x0=x;endtable=[tmp_table];save-asciist_st_eig.dattable;if(~strcmp(getenv('OMIT_PLOTS'),'true'))% PLOTTINGplot(table(:,2),[table(:,4),table(:,6)]);axis([280,460,-1,1]);figure()plot(table(:,3),[table(:,4),table(:,6)]);axis([0,1,-1,1]);% TITLEend% PLOTTING