Nonviscously damped vibrating systems are characterized by dissipative mechanisms depending on the time-history of the response velocity, introduced in the physical models using convolution integrals involving hereditary kernel functions. One of the most used damping viscoelastic models is the Biot's model, whose hereditary functions are assumed to be exponential kernels. The free-motion equations of these types of nonviscous systems lead to a nonlinear eigenvalue problem enclosing certain number of the so-called nonviscous modes with nonoscillatory nature. Traditionally, the nonviscous modes (eigenvalues and eigenvectors) for nonproportional systems have been computed using the state-space approach, computationally expensive. In this paper, we address this problem developing a new method, computationally more efficient than that based on the state-space approach. It will be shown that real eigenvalues and eigenvectors of viscoelastically damped system can be obtained from a linear eigenvalue problem with the same size as the physical system. The numerical approach can even be enhanced to solve highly damped problems. The theoretical results are validated using a numerical example.