Unified constitutive equations have been developed to model the behaviour of metallic materials under various processing conditions. These constitutive equations usually take the form of a set of ordinary differential equations (ODEs), which must be solved thousands of times in a finite element (FE) process simulation. Thus, an efficient and reliable numerical integration method for large systems is crucial for solving this problem. However, in many constitutive equations, numerical stiffness is often present. This means that the stability requirements, rather than the accuracy, constrain the step size. Therefore, certain numerical methods become unsuitable when the required step size becomes unacceptably small. In this study, a series of mathematical analyses was performed to investigate the difficulties in the numerical integration of three sets of unified viscoplastic/creep constitutive equations. Based on an analysis of the current stiffness assessment methods, a novel index was introduced, that enables an accurate assessment of the stiffness of the ODE-type unified constitutive equations. A computational study was also conducted to benchmark several promising implicit numerical integration methods for viscoplastic/creep constitutive equations. This study can assist researchers in metal forming and other fields in choosing appropriate numerical methods when dealing with stiff ODEs.