A 3D multidomain pseudospectral time-domain method is developed for elastic wave equations. The method is based on the spectral derivative operator approximated by Chebyshev or Lagrange polynomials. Unlike the Fourier method that assumes periodic boundary conditions, the Chebyshev pseudospectral method allows for the incorporation of various boundary conditions (such as the free surface boundary condition) into the numerical scheme. In this multidomain scheme, the computational domain is decomposed into a set of subdomains conformal to the problem geometry. Each curved subdomain is then mapped onto a cube in the curvilinear coordinates so that a tensor-product Chebyshev grid can be utilized without the staircasing error. An unsplit-field perfectly matched layer is developed as the absorbing boundary condition. Numerical examples show that this scheme is efficient for simulating elastic waves phenomena in the presence of complex objects. The method is found to be significantly more efficient than the finite-difference time-domain method in terms of memory and run-time requirements.