% example of integrating two coupled ODE's with matlab using Euler's
% method
% A -> B -> C irreversible rxn in isothermal, isobaric PFR

k1 = 1; % rate constant of A -> B 
k2 = 1; % rate constant of B -> C 

% specify initial conditions
CA(1) = 1;
CB(1) = 0;
tau(1) = 0;

% specify step size
dtau = 0.1;

for i = 1:24
    dCAdt = -k1*CA(i);
    dCBdt = k1*CA(i) - k2*CB(i);
    CA(i+1) = CA(i) + dCAdt*dtau;
    CB(i+1) = CB(i) + dCBdt*dtau;
    tau(i+1) = tau(i) + dtau;
end

% now do post-integration calculations
Aprice = 1;
Bprice = 1;
cost = Aprice - Bprice*CB; 
% BE SURE TO USE DOT OPERATIONS HERE IF NEED TO MULTIPLY OR DIVIDE
% TWO VECTORS BY INDIVIDUAL ELEMENT, E.G., CA .* CB 

plot(tau,CA,'b',tau,CB,'r',tau,cost,'k--')
title('conc vs. tau, blue = CA, red = CB, dashed = cost')
ylabel('conc')
xlabel('t')