A simple computer code is presented for simulating the dependence of sulfate reduction on sulfate concentration using Monod kinetics. Unlike previous models, the code provides a numerical initial value problem solution, rather than a two-point boundary value solution, for the Monod model using a search procedure to find the correct starting value for the derivative of sulfate concentration with respect to depth. Accordingly the code is not restricted to cases where sulfate vanishes at finite depth but also can be used to model situations where organic matter is exhausted before total depletion of sulfate can occur. In such situations, the code demonstrates that profiles generated using Monod kinetics differ significantly from those generated using the simple sulfate-independent model proposed by Berner (1964), even when the asymptotic concentration of sulfate at depth remains well above the Monod saturation constant.