In order to model realistic quantum devices it is necessary to simulate quantum systems strongly coupled to their environment. To date, most understanding of open quantum systems is restricted either to weak system–bath couplings or to special cases where specific numerical techniques become effective. Here we present a general and yet exact numerical approach that efficiently describes the time evolution of a quantum system coupled to a non-Markovian harmonic environment. Our method relies on expressing the system state and its propagator as a matrix product state and operator, respectively, and using a singular value decomposition to compress the description of the state as time evolves. We demonstrate the power and flexibility of our approach by numerically identifying the localisation transition of the Ohmic spin-boson model, and considering a model with widely separated environmental timescales arising for a pair of spins embedded in a common environment. 