Ultrasound-based shear wave elastography is an emerging diagnostic procedure that provides estimates of tissue stiffness. Externally induced time-harmonic shear waves enable deeper imaging, enhancing diagnostic accuracy with complete-organ elastograms. Most existing processing techniques rely on time-of-flight, local frequency estimates, and phase-gradient methods. However, these approaches struggle for multidirectional wavefields and require preprocessing steps, which may reduce resolution. Our strategy instead leverages the physics of wave propagation through a Full-Waveform Inversion (FWI) workflow. We demonstrate that FWI has great potential in time-harmonic elastography (THE). It iteratively minimizes an objective function based on the discrepancy between measurements and synthetics. We propose several numerical strategies, such as accumulating multi-frequency measurements and sequentially combining inversion results, rather than incorporating all frequency data into a single objective functional. We show for both recorded and synthetic data, the resulting FWI-based shear-wave speed estimates have higher accuracy, outperforming those from the phase-gradient method.