Viscoelasticity has been widely studied over the past few decades as a combination of the viscous effects of energy dissipation and the elastic effects of energy storage. As an application, viscoelastic dampers are commonly used to suppress the dynamic behaviour of structures. However, in an actual physical environment, it is difficult to explore the stationary probability densities of random responses for viscoelastic systems. For this purpose, we present a numerical method to investigate the dynamic behaviour of a generalized Maxwell-type viscoelastic system under harmonic and Gaussian white noise excitations. Using approximate equivalent and stochastic averaging, we establish an averaged Itô differential equation for the amplitude of the system. For primary external resonance, we solve the reduced Fokker-Planck-Kolmogorov equation by using the successive over-relaxation technique combined with a finite difference method. Through Monte Carlo simulations, we verify the applicability, accuracy and efficiency of the proposed methodology and demonstrate that viscoelasticity has a significant impact on the dynamic behaviour of the viscoelastic systems. This work reveals the remarkable influences of viscoelasticity, excitation intensities, linear and nonlinear damping and linear and nonlinear stiffness on the probability density functions of system responses, which can guide the selection of materials, stiffnesses and structures in viscoelastic damper design.