This article treats the computation of discretized constitutive models of evolutionary-type (like models of viscoelasticity, plasticity, and viscoplasticity) with quasi-static finite elements using diagonally implicit Runge-Kutta methods (DIRK) combined with the Multilevel-Newton algorithm (MLNA). The main emphasis is on promoting iterative methods, as opposed to the more traditional direct methods, for solving the non-symmetric systems which occur within the DIRK/MLNA. It is shown that iterative solution of the arising sequences of linear systems can be substantially accelerated by various techniques that aim at sharing part of the computational effort throughout the sequence. In this way iterative solution becomes attractive at clearly lower dimensions than the dimensions where direct solvers start to fail for memory reasons. The applications are related to small strain viscoplasticity of a smooth constitutive model for plastics and a finite strain plasticity model with non-linear kinematic hardening developed for metals.