The inversion of teleseismic receiver functions for lithospheric structure is difficult due to the non-linearity of the problem, which is greatly increased in the presence of dipping interfaces and layer anisotropy. Given an efficient ray-theoretical tool for forward-modelling teleseismic seismograms, we perform a directed Monte Carlo search technique using the neighbourhood algorithm of Sambridge, enabling us to search 20–30 parameters in a reasonable amount of computer time. Tests on synthetic data reveal inherent velocity–depth trade-offs in typical data sets, due to the limited moveout present in teleseismic Ps; the azimuth of the anisotropic symmetry axis and the strike of a dipping interface prove to be well-resolved given adequate backazimuthal coverage. We apply this technique to two single-station data sets. The first, from permanent station PGC, Vancouver Island, British Columbia, displays dipping low-velocity sediment layers in the mid-crust. The second, from a station at the northern end of the Tibetan plateau operating in 1991 and 1992, requires a sequence of thick crustal anisotropic layers to explain the observed pattern of receiver-function arrivals.