We introduce a matrix-product state based method to efficiently obtain dynamical response functions for two-dimensional microscopic Hamiltonians, which we apply to different phases of the Kitaev-Heisenberg model. We find significant broad high energy features beyond spin-wave theory even in the ordered phases proximate to spin liquids. This includes the phase with zig-zag order of the type observed in α-RuCl3, where we find high energy features like those seen in inelastic neutron scattering experiments. Our results provide an example of a natural path for proximate spin liquid features to arise at high energies above a conventionally ordered state, as the diffuse remnants of spin-wave bands intersect to yield a broad peak at the Brillouin zone center.Introduction. The interplay of strong interactions and quantum fluctuations in spin systems can give rise to new and exciting physics. A prominent example are quantum spin liquids (QSL), as fascinating as they are hard to detect: they lack local order parameters and are instead characterized in terms of emergent gauge fields. On the experimental side, spectroscopic measurements provide particularly useful insights into such systems, in particular by probing the fractionalised excitations (e.g. deconfined spinons) accompanying the gauge field. Such measurements can be related to dynamical response functions, e.g. inelastic neutron scattering to the dynamical structure factor. On the theoretical side, determining the ground state properties of such quantum spin models is already a hard problem, and it is even more challenging to understand the dynamics of local excitations.Here we present a combination of the density-matrix renormalization (DMRG) ground state method and a matrix-product states (MPS) based dynamical algorithm to obtain the response functions for generic twodimensional spin systems. With this we are able to access the dynamics of exotic phases that can occur in frustrated systems. Moreover it is also very useful for regular ordered phases where one would conventionally use large-S approximations, which in some cases cannot qualitatively explain certain high energy features 1,2 . We demonstrate our method by applying it to the currently much-studied Kitaev-Heisenberg model (KHM) model on the honeycomb lattice