The compressible Navier-Stokes equations can be extended to model multi-species, chemically reacting gas flows. The result is a large system of convection-diffusion equations with stiff source terms. In this paper, we develop the framework needed to apply modern high accuracy numerical methods from computational gas dynamics to this extended system. We implement these developments with a particular high accuracy characteristic-based method, the finite difference ENO space discretization with the third-order TVD Runge-Kutta time discretization, combined with the second-order accurate Strang time splitting of the reaction terms. We illustrate the capabilities of this approach with calculations of a one-dimensional reacting shock tube and a two-dimensional combustor.