To decouple the parameters in elastic full-waveform inversion (FWI), we evaluated a new multistep-length gradient approach to assign individual weights separately for each parameter gradient and search for an optimal step length along the composite gradient direction. To perform wavefield extrapolations for the inversion, we used parallelized high-precision finite-element (FE) modeling in the time domain. The inversion was implemented in the frequency domain; the data were obtained at every subsurface grid point using the discrete Fourier transform at each time-domain extrapolation step. We also used frequency selection to reduce cycle skipping, time windowing to remove the artifacts associated with different source spatial patterns between the test and predicted data, and source wavelet estimation at the receivers over the full frequency spectrum by using a fast Fourier transform. In the inversion, the velocity and density reconstructions behaved differently; as a low-wavenumber tomography (for velocities) and as a high-wavenumber migration (for density). Because velocities and density were coupled to some extent, variations were usually underestimated (smoothed) for and and correspondingly overestimated (sharpened) for . The impedances and from the products of the velocity and density results compensated for the under- or overestimations of their variations, so the recovered impedances were closer to the correct ones than , , and were separately. Simultaneous reconstruction of , , and was robust on the FE and finite-difference synthetic data (without surface waves) from the elastic Marmousi-2 model; satisfactory results are obtained for , , , and the recovered and from their products. Convergence is fast, needing only a few tens of iterations, rather than a few hundreds of iterations that are typical in most other elastic FWI algorithms.