function yprime = chemreac(t,y);
% chemistry dynamics posed in spreadsheet book
% return a column vector
% note - the variables k1, km1 were declare 'global' in the driver
% y(1)=CO; y(2)=H2O; y(3)=CO2; y(4)=H2
% since CO and H2O satisfy same differential equation, only compute
% term once
yp1 = -k1*y(1)*y(2) + km1*y(3)*y(4);
% same for CO2 and H2
yp3 = -km1*y(3)*y(4) + k1*y(1)*y(2);
yprime = [yp1; yp1; yp3; yp3];