% solve differential equation
% for isothermal batch reactor
% using ode15s as a DAE system
k = 14.7 ; % 1/(mol/L)/min
Ca0 = 0.0209 ; % mol/L
Cb0 = Ca0/3 ;
Cc0 = 0
tspan = 0:0.1:20 ; % solution times
% set mass matrix
M = [ 1 0 0 ; % 1 on diagonal identifies ODE
      0 0 0 ; % 0 on diagonal identifies algebraic
      0 0 0 ];
y0 = [ Ca0 ; Cb0 ; Cc0 ]; % initial conditions
options = odeset('Mass',M); % set for DAE
% anonymous function to include extra arguments
batchanon = @(t,y) batchDAE(t,y,k,Ca0,Cb0);
% solve system with ode15s
[ t, ysoln ] = ode15s(batchanon,tspan,y0,options);
% unpack solution
Ca = ysoln(:,1); Cb = ysoln(:,2) ; Cc = ysoln(:,3);
% create plot
plot(t,Ca,'b',t,Cb,'m',t,Cc,'g')
grid
xlabel('Time - minutes')
ylabel('Concentration - mol/L')
title('Batch Reactor DAE Solution')
legend('Ca','Cb','Cc')

