% Multiple reactions simulation
% using ode15s and a DAE format
% rate constants
k1 = 14.7; % 1/(mol/L)/min
k2 = 1.53;
k3 = 0.294;
% initial conditions
A0 = 0.0209; % mol/L
B0 = A0/3;
C0 = 0;
D0 = 0;
E0 = 0;
F0 = 0;
y0 = [ A0 B0 C0 D0 E0 F0 ]';
% time span
tspan = 0:500;
% mass matrix
M = [ 1 0 0 0 0 0 ;
      0 1 0 0 0 0 ;
      0 0 1 0 0 0 ;
      0 0 0 1 0 0 ;
      0 0 0 0 0 0 ;
      0 0 0 0 0 0 ];
% set options for DAE
options = odeset('Mass',M);
% anonymous function for the DAE system
multibatchanon = @(t,y) multibatch(t,y,k1,k2,k3,A0);
% solve system
[ t,y ] = ode15s(multibatchanon,tspan,y0,options);
% unpack solution
A = y(:,1);
B = y(:,2);
C = y(:,3);
D = y(:,4);
E = y(:,5);
F = y(:,6);
% create plot
plot(t,A,'b-',t,B,'r-',t,C,'m-',t,D,'k-',t,E,'c-',t,F,'g-');
grid
axis([ 0 500 0 0.025 ]);
legend('A','B','C','D','E','F','Location','northeastoutside')
xlabel('Time - minutes')
ylabel('Concentration - mol/L')
title('Isothermal Batch Reactor - Multiple Reactions')
