The one-dimensional contact process (CP) in a heterogeneous environment-a binary chain consisting of two types of site with different recovery rates-is investigated. It is argued that the commonly used random-sequential Monte Carlo simulation method which employs a discrete notion of time is not faithful to the rates of the contact process in a heterogeneous environment. Therefore, a modification of this algorithm along with two alternative continuous-time implementations are analyzed. The latter two are an adapted version of the n -fold way used in Ising model simulations and a method based on a modified priority queue. It is demonstrated that the commonly used (but incorrect as we believe) discrete-time method yields a different critical threshold from all other algorithms considered. Finite-size scaling of the lowest gap in the spectrum of the Liouville time-evolution operator for the CP gives an estimate of the critical rate which supports these findings. Further, a performance test indicates an advantage in using the continuous-time methods in systems with heterogeneous rates. This result promises to help in the analysis of the CP in disordered systems with heterogeneous rates in which simulation is a challenging task due to very long relaxation times.