In the absence of a full analytical treatment of nonlinear structure formation in the universe, numerical simulations provide the critical link between the properties of the underlying model and the features of the observed structures. Currently N-body simulations are the main tool to study structure growth. We explore an alternative framework for numerical simulations of structure formation. The underlying idea is to replace the long-range gravitational force in the Vlasov-Poisson system by a purely local interaction. To this end we trade the classical phase space distribution for its quantum mechanical counterpart, the Wigner distribution function. Its dynamical equation is equivalent to the Schrödinger equation and reduces to the Vlasov equation in the formal semi-classical limit. The proposed framework allows in principle to simulate systems with arbitrary phase space distributions and could for instance be beneficial for simulations of warm dark matter, where the velocity dispersion is important. We discuss several methods to obtain a set of wavefunctions whose Wigner distribution is close to a given initial phase space distribution function. An auxiliary gauge field is introduced to mediate the gravitational interaction, thereby obtaining a local Schrödinger-Maxwell system. We also use the ideas of lattice gauge theories to obtain a fully gauge-invariant discretization of the equations of motion. Their iterative solution was implemented in a three-dimensional simulation code. We discuss its computational complexity and memory requirements. Several testbed simulations were performed with this method. We compared the gravitational collapse of a Gaussian wavefunction with an independent numerical solution of the spherically symmetric Schrödinger-Newton system. The results were found to be in good agreement. Finally, a simple example of the growth of cosmic perturbations is investigated within our framework. We conclude by outlining various possible directions to optimize and develop our method.