We investigate the ability of basis function expansions to reproduce the evolution of a Milky Way-like dark matter halo, extracted from a cosmological zoom-in simulation. For each snapshot, the density of the halo is reduced to a basis function expansion, with interpolation used to recreate the evolution between snapshots. The angular variation of the halo density is described by spherical harmonics, and the radial variation either by biorthonormal basis functions adapted to handle truncated haloes or by splines. High fidelity orbit reconstructions are attainable using either method with similar computational expense. We quantify how the error in the reconstructed orbits varies with expansion order and snapshot spacing. Despite the many possible biorthonormal expansions, it is hard to beat a conventional Hernquist-Ostriker expansion with a moderate number of terms ($gtrsim15$ radial and $gtrsim6$ angular). As two applications of the developed machinery, we assess the impact of the time-dependence of the potential on (i) the orbits of Milky Way satellites, and (ii) planes of satellites as observed in the Milky Way and other nearby galaxies. Time evolution over the last 5 Gyr introduces an uncertainty in the Milky Way satellites orbital parameters of $sim 15$ per cent, comparable to that induced by the observational errors or the uncertainty in the present-day Milky Way potential. On average, planes of satellites grow at similar rates in evolving and time-independent potentials. There can be more, or less, growth in the planes thickness, if the plane becomes less, or more, aligned with the major or minor axis of the evolving halo.