[en] Highlights: • A new method for solving chemical reaction kinetics with disorder in the rates (unimolecular or bimolecular) is presented. • The method uses a mean-field approximation to quickly compute self-consistent steady state populations. • We illustrate the proposed method with a working example of H2 production on a heterogeneous surface. We propose a method to quickly compute steady state populations of species undergoing a set of chemical reactions whose rate constants are heterogeneous. Using an average environment in place of an explicit nearest neighbor configuration, we obtain a set of equations describing a single fluctuating active site in the presence of an averaged bath. We apply this Mean Field Steady State (MFSS) method to a model of $H_2$ production on a disordered surface for which the activation energy for the reaction varies from site to site. The MFSS populations quantitatively reproduce the KMC results across the range of rate parameters considered.