% cstr_run.m
% b.w. bequette
% 2 Dec 96
% revised - 4 Dec 96 - more documentation
%
% generate phase-plane plot for nonlinear cstr (cstr.m)
%
% set of initial conditions
%
  x01a= [0.5;305];
  x01 = [0.5;315];
  x02 = [0.5;325];
  x03 = [0.5;335];
  x04 = [0.5;345];
  x05 = [0.5;355];
  x06 = [0.5;365];
  x07 = [0.5;375];
  x08 = [0.5;385];
  x09 = [0.5;395];
  x010= [0.5;405];
  x011= [0.5;415];
  x012a=[9.5;305];
  x012= [9.5;315];
  x013= [9.5;325];
  x014 = [9.5;335];
  x015 = [9.5;345];
  x016 = [9.5;355];
  x017 = [9.5;365];
  x018 = [9.5;375];
  x019 = [9.5;385];
  x020 = [9.5;395];
  x021 = [9.5;405];
  x022 = [9.5;415];
%
% use ode45 for numerical integration
%
  [t1a,x1a] = ode45('cstr',0,10,x01a);
  [t1,x1] = ode45('cstr',0,10,x01);
  [t2,x2] = ode45('cstr',0,10,x02);
  [t3,x3] = ode45('cstr',0,10,x03);
  [t4,x4] = ode45('cstr',0,10,x04);
  [t5,x5] = ode45('cstr',0,10,x05);
  [t6,x6] = ode45('cstr',0,10,x06);
  [t7,x7] = ode45('cstr',0,10,x07);
  [t8,x8] = ode45('cstr',0,10,x08);
  [t9,x9] = ode45('cstr',0,10,x09);
  [t10,x10] = ode45('cstr',0,10,x010);
  [t11,x11] = ode45('cstr',0,10,x011);
  [t12,x12] = ode45('cstr',0,10,x012);
  [t12a,x12a] = ode45('cstr',0,10,x012a);  
  [t13,x13] = ode45('cstr',0,10,x013);
  [t14,x14] = ode45('cstr',0,10,x014);
  [t15,x15] = ode45('cstr',0,10,x015);
  [t16,x16] = ode45('cstr',0,10,x016);
  [t17,x17] = ode45('cstr',0,10,x017);
  [t18,x18] = ode45('cstr',0,10,x018);
  [t19,x19] = ode45('cstr',0,10,x019);
  [t20,x20] = ode45('cstr',0,10,x020);
  [t21,x21] = ode45('cstr',0,10,x021);
  [t22,x22] = ode45('cstr',0,10,x022);
%
  x040 = [0.5;370];
  x041 = [0.5;372.5];
  x049 = [9.5;320];
  x050 = [9.5;326.25];
  x051 = [9.5;327.5];
  x052 = [9.5;330];
  x053 = [9.5;332.5];
  [t40,x40] = ode45('cstr',0,10,x040);
  [t41,x41] = ode45('cstr',0,10,x041);
  [t49,x49] = ode45('cstr',0,10,x049);
  [t50,x50] = ode45('cstr',0,10,x050);
  [t51,x51] = ode45('cstr',0,10,x051);
  [t52,x52] = ode45('cstr',0,10,x052);
  [t53,x53] = ode45('cstr',0,10,x053);
%
% phase-plane plot
%
  plot(x1(:,1),x1(:,2),'w')
  hold on
  plot(x1a(:,1),x1a(:,2),'w')
  plot(x2(:,1),x2(:,2),'w')
  plot(x3(:,1),x3(:,2),'w')
  plot(x4(:,1),x4(:,2),'w')
  plot(x5(:,1),x5(:,2),'w')
  plot(x6(:,1),x6(:,2),'w')
  plot(x7(:,1),x7(:,2),'w')
  plot(x8(:,1),x8(:,2),'w')
  plot(x9(:,1),x9(:,2),'w')
  plot(x10(:,1),x10(:,2),'w')
  plot(x11(:,1),x11(:,2),'w')
  plot(x12(:,1),x12(:,2),'w')
  plot(x12a(:,1),x12a(:,2),'w')
  plot(x13(:,1),x13(:,2),'w')
  plot(x14(:,1),x14(:,2),'w')
  plot(x15(:,1),x15(:,2),'w')
  plot(x16(:,1),x16(:,2),'w')
  plot(x17(:,1),x17(:,2),'w')
  plot(x18(:,1),x18(:,2),'w')
  plot(x19(:,1),x19(:,2),'w')
  plot(x20(:,1),x20(:,2),'w')
  plot(x21(:,1),x21(:,2),'w')
  plot(x22(:,1),x22(:,2),'w')
% 
% mark the stable and unstable steady-states (o = stable, + = unstable)
%
  plot([8.564],[311.2],'wo',[5.518],[339.1],'w+',[2.359],[368.1],'wo')
%
%  some initial conditions not used
%  x030 = [5.768;339.1];
%  x031 = [5.368;339.1];
%  x032 = [5.518;338.5];
%  x033 = [5.518;339.7];
%  [t30,x30] = ode45('cstr',0,10,x030);
%  [t31,x31] = ode45('cstr',0,10,x031);
%  [t32,x32] = ode45('cstr',0,10,x032);
%  [t33,x33] = ode45('cstr',0,10,x033);
%
  plot(x40(:,1),x40(:,2),'w');
  plot(x41(:,1),x41(:,2),'w');
  plot(x49(:,1),x49(:,2),'w');
  plot(x50(:,1),x50(:,2),'w');
  plot(x51(:,1),x51(:,2),'w');
  plot(x52(:,1),x52(:,2),'w');
  plot(x53(:,1),x53(:,2),'w');
%
% plot(x30(:,1),x30(:,2),'w');
% plot(x31(:,1),x31(:,2),'w');
% plot(x32(:,1),x32(:,2),'w');
% plot(x33(:,1),x33(:,2),'w');
%
% set limits on plot, and label axes and provide plot title
%
  axis([0 10 300 425])
  xlabel('concentration, kgmol/m^3')
  ylabel('temperature, deg K')
  title('cstr with multiple steady-states')