The aim of this work is to propose a mathematical model in terms of an exact analytical solution that may be used in numerical simulation and prediction of oscillatory dynamics of a one-dimensional viscoelastic system experiencing large deformations response. The model is represented with the use of a mechanical oscillator consisting of an inertial body attached to a nonlinear viscoelastic spring. As a result, a second-order first-degree Painlevé equation has been obtained as a law, governing the nonlinear oscillatory dynamics of the viscoelastic system. Analytical resolution of the evolution equation predicts the existence of three solutions and hence three damping modes of free vibration well known in dynamics of viscoelastically damped oscillating systems. Following the specific values of damping strength, over-damped, critically-damped and under-damped solutions have been obtained. It is observed that the rate of decay is not only governed by the damping degree but, also by the magnitude of the stiffness nonlinearity controlling parameter. Computational simulations demonstrated that numerical solutions match analytical results very well. It is found that the developed mathematical model includes a nonlinear extension of the classical damped linear harmonic oscillator and incorporates the Lambert nonlinear oscillatory equation with well-known solutions as special case. Finally, the three damped responses of the current mathematical model devoted for representing mechanical systems undergoing large deformations and viscoelastic behavior are found to be asymptotically stable.