Seismic anisotropy, wave attenuation, and dispersion are critical phenomena of wave propagation in real media. Full-wavefield modeling of wave behavior in such media plays an important role in investigating dynamic features of the earth’s interior and in full-waveform inversion of anisotropic parameters and velocity dispersion. We have developed a numerical scheme to model the full-waveform response from a point source in a vertically varying viscoelastic medium of arbitrary anisotropy. The method is implemented in the frequency domain, so the complexity of anelasticity and anisotropy can be described by a complex elastic stiffness matrix and frequency-dependent moduli can also be readily incorporated. In our scheme, we solve the elastodynamic equations for general anisotropy through finite-element method in the frequency-wavenumber domain and we use the stiffness reduction method to suppress reflections from artificial boundaries along the depth direction. A nonuniform 2D Fourier transform strategy is developed to reconstruct the spatial-domain counterparts from the wavenumber-domain solutions. The time-domain responses are then obtained by taking inverse fast Fourier transform with respect to frequency. We validate the method by comparing the numerical results with the exact solutions for a homogeneous transversely isotropic model and a two-layered model. In the application example, we further proved the feasibility and generality of the scheme using an attenuative, dispersive model with velocity and attenuation anisotropy. Our scheme enjoys significant advantages in incorporating various viscoelastic/dispersive behaviors and general anisotropy, and thus provides a useful tool for numerical simulation of dynamic responses in practical application.