On the numerical solution of a BGK-type model for chemical reactions