It is now widely accepted that elastic properties of the continental lithosphere and the underlying sublithospheric mantle are both anisotropic and laterally heterogeneous at a range of scales. To fully exploit modern three-component broad-band array data sets requires the use of comprehensive modelling tools. In this work, we investigate the use of a wide-angle, one-way wave equation to model variations in teleseismic 3-D waveforms due to 2-D elastic heterogeneity and anisotropy. The one-way operators are derived based on a high-frequency approximation of the square-root operator and include the effects of wave propagation as well as multiple scattering. Computational cost is reduced through a number of physically motivated approximations. We present synthetic results from simple 1-D (layer over a half-space) and 2-D (subduction zone) models that are compared with reference solutions. The algorithm is then used to model data from an array of broad-band seismograph stations deployed in northwestern Canada as part of the IRIS-PASSCAL/LITHOPROBE CANOE experiment. In this region radial-component receiver functions show a clear continental Moho and the presence of crustal material dipping into the mantle at the suture of two Palaeo-Proterozoic terranes. The geometry of the suture is better defined on the transverse component where subduction is associated with a ∼10 km thick layer exhibiting strong elastic anisotropy. The modelling reproduces the main features of the receiver functions, including the effects of anisotropy, heterogeneity and finite-frequency scattering.