An unconditionally stable numerical procedure is developed to integrate the stiff differential equation associated with the flow rule for viscoplasticity used by Swegle and Grady [J. Appl. Phys.58, 692 (1985)] to predict shock structure. Following the work of Rubin [J. Appl. Math. Phys. ZAMP40, 846 (1989)] the flow rule is integrated implicitly and is converted into a scalar equation. It is shown that for the Swegle–Grady model this scalar equation reduces to a quadratic equation that can be solved without iteration.