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

c0 = [0 333.15]';
[t, c] = ode15s('model43', [0 2000], 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;
clear c t C T;

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


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

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

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

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

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

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

c0 = [8 273.15]';
[t,c] = ode15s('model43', [0 200], c0);
C = c(:,1);
T=c(:,2) -273.15;
plot (T,C);
hold on;
clear c t C T;

c0 = [14 273.15]';
[t,c] = ode15s('model43', [0 200],c0);
C = c(:,1);
T=c(:,2) -273.15;
plot (T,C);
hold on;
clear c t C T;

c0 = [14 323.15]';
[t,c] = ode15s('model43', [0 200],c0);
C = c(:,1);  T = c(:,2) - 273.15;
plot(T,C);
hold on;
clear c t C T;

c0 = [10 323.15]';
[t,c] = ode15s('model43', [0 200],c0);
C = c(:,1);  T = c(:,2) - 273.15;
plot(T,C);
hold on;
clear c t C T;

c0 = [8 323.15]';
[t,c] = ode15s('model43', [0 200], c0);
C = c(:,1);  T = c(:,2) - 273.15;
plot(T,C);
