One of the successful approaches to model the time-dependent behaviour of elastomers is proposed by Bergstrom and Boyce (JMPS 46:931-954, 1998). The model is micromechanically inspired from the relaxation of a single entangled chain in a polymer gel matrix. Although the theory of inelasticity based on multiplicative decomposition of the deformation gradient is well established, the complexity of the nonlinear evolution law as well as the nonlinear equilibrium and non-equilibrium material response necessitates a precise description of the algorithmic setting. This contribution presents for the first time a novel numerical implementation of the Bergstrom-Boyce model in the context of finite element analysis and elaborates theoretical aspects of the model. The thermodynamical consistency of the evolution law is proven and a parameter study with respect to the material parameters has been carried out. The agreement of the model with the recent experimental data is investigated.