ABSTRACT We describe and test a family of new numerical methods to solve the Schrödinger equation in self-gravitating systems, e.g. Bose–Einstein condensates or ‘fuzzy’/ultra-light scalar field dark matter. The methods are finite-volume Godunov schemes with stable, higher order accurate gradient estimation, based on a generalization of recent mesh-free finite-mass Godunov methods. They couple easily to particle-based N-body gravity solvers (with or without other fluids, e.g. baryons), are numerically stable, and computationally efficient. Different sub-methods allow for manifest conservation of mass, momentum, and energy. We consider a variety of test problems and demonstrate that these can accurately recover solutions and remain stable even in noisy, poorly resolved systems, with dramatically reduced noise compared to some other proposed implementations (though certain types of discontinuities remain challenging). This is non-trivial because the ‘quantum pressure’ is neither isotropic nor positive definite and depends on higher order gradients of the density field. We implement and test the method in the code gizmo.