function dy = chem_kinetic_eqn(t,y)
global k1 km1 k2 
dy = zeros(4,1); % need to be a column vector
dy(1) = -k1*y(1)*y(2)+km1*y(3);
dy(2) = -k1*y(1)*y(2)+km1*y(3)+k2*y(3);
dy(3) = k1*y(1)*y(2)-km1*y(3)-k2*y(3);
dy(4) = k2*y(3);