We present a method for modelling seismic wave propagation in a whole-earth model by solving the elastodynamic equations in 2-D cylindrical coordinates (r, 0) using the Fourier pseudospectral method (PSM). In solving the 2-D cylindrical elastodynamic equations for a whole-earth model, a singularity arises at the centre (r = 0) of the earth. To avoid the singularity, we develop a scheme that uses extension of field variables in the radial direction, with which computation of the wavefield at the centre is avoided, so that the wave propagation through the centre can be calculated. The time interval used in the calculation is determined by the smallest lateral grid spacing around the centre in the model. In a cylindrical coordinate system, the smallest lateral grid spacing is generally so small that the calculation is too time-consuming to be realistically carried out even on a supercomputer. We adopt a multidomain scheme to increase the smallest lateral grid spacing and avoid the oversampling of the physical domain around the centre of the earth. A smoothing scheme in the wavenumber domain is also proposed, which enables us to use a large enough time interval to allow the calculation for the whole-earth model on a desktop workstation. The waveforms calculated by the present method are compared with those obtained by the Direct Solution Method (DSM) to demonstrate their high accuracy. This method significantly reduces the computer memory and computation time required and makes it possible to study the effects of small-wavelength heterogeneities that can be approximated as azimuthally symmetric on wave propagation in the earth. We apply the present method to study the effects of local heterogeneity in the earth by adding a low-velocity perturbation above the core-mantle boundary (CMB) to the IASP91 earth model.