Ultrasonic imaging system evaluation is often based on models of the transducer as a distribution of baffled piston sources, and of the tissue as a homogeneous, linear acoustic medium, e.g., Jensen’s Field code. In reality, these are fairly gross idealizations, since the transducer exhibits more complicated response modes and real tissue is inhomogeneous and nonlinear. Greater model fidelity would be useful, especially in the context of transducer design qualification, second harmonic imaging, and acoustic power indices. To this end we combine 2D finite element models of transducer dynamics with highly accurate 2D finite difference propagation models in the large-scale inhomogeneous tissue cross-sections. Transducer models employ the time-domain code, PZFlex, and tissue models utilize a new pseudospectral solver to be included in PZFlex. The pseudospectral algorithm solves the inhomogeneous acoustic wave equation using FFTs for high order approximation of the spatial differential operator and a fourth-order, explicit time integrator. Second-order (B/A) nonlinearity and frequency-accurate, causal absorption are included. We describe the algorithmic and modeling issues, and present a suite of simulations in lossy, nonlinear abdominal cross sections and tissue showing coupling of the 1D medical array to the tissue model and scattering from deeper inhomogeneities and back to the transducer. In contrast to paraxial schemes, like the KZK method, details of the field transmitted from the transducer and all backscatter within the model are included. However, models are currently limited to 2D (plane or axisymmetric) on readily available hardware.