Simulation of Rayleigh waves requires high accuracy and an adequate spatial sampling at the surface. Discrete cosine and sine transforms are used to compute spatial derivatives along the direction perpendicular to the surface of the earth. Unlike the standard Fourier method, these transforms allow nonperiodic boundary conditions to be satisfied, in particular, the stress-free conditions at the surface. Because simulation of surface waves requires more points per minimum wavelength at the surface than simulation of body waves, the equispaced grid is not efficient. To overcome this problem, a grid compression is performed at the surface to obtain a denser spatial sampling. Grid size is minimal at the surface and increases with depth until reaching, at a relatively shallow depth, the grid points per wavelength required by the body waves. The stress-free boundary conditions are naturally handled by expanding the appropriate stress components in terms of the discrete sine transform. The wave equation is solved in the particle-velocity and stress formulation using a Runge-Kutta time integration and the convolutional PML (CPML) method to prevent reflections from the mesh boundaries. The simulations are very accurate for shallow sources and receivers and large offsets.