This paper describes an efficient approach for computing the frequency response of seismic waves propagating in 2- and 3-D earth models within which the magnitude and phase are required at many locations. The approach consists of running an explicit finite difference time domain (TD) code with a time harmonic source out to steady-state. The magnitudes and phases at locations in the model are computed using phase sensitive detection (PSD). PSD does not require storage of time-series (unlike a fast Fourier transform), reducing its memory requirements. Additionally, the response from multiple sources can be obtained from a single finite difference run by encoding each source with a different frequency. For 2-D models with many sources, this time domain phase sensitive detection (TD-PSD) approach has a higher arithmetic complexity than direct solution of the finite difference frequency domain (FD) equations using nested dissection re-ordering (FD-ND). The storage requirements for 2-D finite difference TD-PSD are lower than FD-ND. For 3-D finite difference models, TD-PSD has significantly lower arithmetic complexity and storage requirements than FD-ND, and therefore, may prove useful for computing the frequency response of large 3-D earth models.