% This file calls the stiff differential equation solver ode23s, 
% to solve the autocatalytic batch reactor problem in section 4.2;
%

y0 = [1.0 0 0 0]';
[T, C] = ode23s('model42', [0 0.66], y0);

%
% convert to dimensional cases
%
k1 = 10^(-4);
t = T/k1/3600;

%
% assume CA0=1(mol/l), therefore the
% concentration values can be treated as dimensional directly.
%

CA = C(:, 1);
CB = C(:, 2);
CC = C(:, 3);
CAB = C(:, 4);

plot(t, CA, '-', t, CB, '--', t, CC, '-.', t, CAB, '.');
axis([0 2 0 1.0]);
title('Autocatalytic Batch Reactor (Example 4.2)');
xlabel('Time (hr)');
ylabel('Concentration (kmol/l)');
text(0.3, 0.8, 'A');
text(0.8, 0.5, 'B');
text(0.65, 0.65, 'C');
text(1.2, 0.15, 'AB');
text(1.4, 0.4, 'alpha=8');
%grid on;
