Elastomeric solid materials typically exhibit a pronounced viscoelastic response. In this paper we consider a large deformation viscoelasticity theory for isotropic elastomeric materials which uses a multi-branch multiplicative decomposition of the deformation gradient. We then describe the numerical implementation of the theory in the open-source finite element program FEniCSx. Several example simulations which demonstrate the capability of the theory and its numerical implementation to model stress-relaxation, creep, stretch-rate sensitivity, hysteresis, damped inertial oscillations, and dynamic column buckling are shown. The source codes for these simulations are provided. The theory and the codes presented in this paper lay the foundation for future extensions of the theory and its numerical implementation to include the effects of coupling with thermal, electrical, and magnetic fields — extensions which are of central importance in modeling the response of soft-active materials.