function dy = zeq(t,y,k1,k2,k3)

% === differential equation =======
dy = zeros(4,1);    % a column vector
dy(1) = -k1*y(1)*y(2) + k2*y(3);
dy(2) = -k1*y(1)*y(2) + k2*y(3) + k3*y(3);
dy(3) = k1*y(1)*y(2) - k2*y(3) - k3*y(3);
dy(4) = k3*y(3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
