% This file calls the built-in ode45, 
% to solve the CSTR problem in section 4.3;
%

c0 = [0 333.15]';
[t, c] = ode23s('model43', [0 200], c0);

C = c(:, 1);
T = c(:, 2) - 273.15;

plot(T, C);
axis([0 80 0 16]);
title('Phase Plane (Example 4.3)');
xlabel('Temperature (degree C)');
ylabel('Concentration (mol/l)');
 text(30, 6, 'K=0.0');
%text(30, 6, 'K=1.0');
hold on;

c=0; t=0; C = 0; T = 0;
%
c0 = [0 273.15]';
[t, c] = ode23s('model43', [0 1480], c0);
C = c(:,1);
T = c(:,2) - 273.15;
plot(T,C);
hold on;
c = 0; t=0; C =0; T=0;

c0 = [2 273.15]';
[t, c] = ode23s('model43', [0 1475], c0);
C = c(:,1); T = c(:,2) - 273.15;
plot(T,C);
hold on;
c = 0; t=0; C=0; T=0;

c0 = [16 303.15]';
[t, c] = ode23s('model43', [0 100], c0);
C = c(:,1); T= c(:,2) - 273.15;
plot(T,C);
hold on;
c = 0; t=0; C=0; T=0;

c0 = [16 313.15]';
[t, c] = ode23s('model43', [0 30], c0);
C = c(:,1); T= c(:,2) - 273.15;
plot(T,C);
hold on;
c =0; t=0; C=0; T=0;

c0 = [0 333.15]';
[t, c] = ode23s('model43', [0 156], c0);
C = c(:,1); T= c(:,2) - 273.15;
plot (T,C);
hold on;
c =0; t=0; C=0; T=0;

c0 = [0 313.15]';
[t, c] = ode23s('model43', [0 187], c0);
C = c(:,1); T= c(:,2) - 273.15;
plot (T,C);
hold on;
c =0; t=0; C=0; T=0;

c0 = [0 343.15]';
[t, c] = ode23s('model43', [0 28], c0);
C = c(:,1); T= c(:,2) -273.15;
plot (T,C);
hold on;