An iterative numerical method is used to simulate the nonequilibrium distributions for nonlinear kinetic systems. This method is much simpler and more stable than the full time-dependent solution of the master equation, yet it yields both the phenomenological rate constants and the nonequilibrium internal-state distributions. The method should be general but it is illustrated here by applying it to the dissociation-recombination problem of a diatomic molecule in excess inert gas. In particular, we studied the dissociation of O2 and the recombination of O atoms in excess Ar at 293, 600, 3000, 4000, and 18000 K. We find that, after an induction period, a quasisteady state is reached for which the phenomenological rate law holds with temperature-dependent but concentration-independent rate coefficients over the entire temperature range from 293 to 18000 K. The phenomenological rate law is valid even when there is appreciable back-reaction, and the phenomenological rate constant can be predicted by an eigenvalue analysis even where nonlinear processes play a significant role over the whole steady kinetics regime. The upper and lower limits of the ratio of forward to reverse phenomenological rates for which the rate coefficients are constant are obtained numerically at each temperature. We also calculate the nonequilibrium internal-state distributions at high temperature, and we find that they can be characterized by an invariant vector in concentration space over the entire range for which the phenomenological rate law applies. This is also shown analytically by applying singular-perturbation techniques. The numerical results also confirm the essential correctness of the predictions of the singular-perturbation method of Brau, Hogarth, and McElwain, even at the very high temperature of 18000 K, where the perturbation parameter is not very small.