<p>The dielectric elastomer actuator can generate nonlinear vibrations under applied voltages. At large deformations, especially near the limiting deformations, the dielectric elastomer shows stiffening phenomenon, which is widely modeled by the Gent model. However, this phenomenon can induce local illness for the corresponding governing ordinary differential equations, and also induce a singularity for the equations. Therefore, the conventional explicit simulation methods may have large simulation errors or even fail as the divergence due to the singularity. Here an implicit Gauss-Legendre collocation method with a variable step size scheme is proposed to study the nonlinear dynamics of two typical dielectric elastomer systems. The scheme of variable step size is based on the consideration of the convergence of iterations and the avoidance of singularity. In the absence of periodic forcing, the method can both become energy preserving and can also avoid the divergence due to the singularity. In the presence of periodic forcing, dynamic behaviors for a large region in the parameter plane cannot be determined due to the divergence when using the classic fourth order Runge-Kutta method. While for the presented Gauss-Legendre collocation method, dynamic behaviors in the whole parameter region are located, and its accuracy is further verified by the global dynamics in the phase plane which include periodic vibrations or chaotic vibrations. Thus, the proposed method will be useful in simulating nonlinear dynamical behavior with stiffening.</p>