The diversity of structures in the Universe (from the smallest galaxies to the largest superclusters) has formed under the pull of gravity from the tiny primordial perturbations that we see imprinted in the cosmic microwave background. A quantitative description of this process would require description of motion of zillions of dark matter particles. This impossible task is usually circumvented by coarse graining the problem: one either considers a Newtonian dynamics of 'particles' with macroscopically large masses or approximates the dark matter distribution with a continuous density field. There is no closed system of equations for the evolution of the matter density field alone and instead it should still be discretized at each time step. In this work, we describe a method of solving the full six-dimensional Vlasov-Poisson equation via a system of auxiliary Schrodinger-like equations. The complexity of the problem gets shifted into the choice of the number and shape of the initial wavefunctions that should only be specified at the beginning of the computation (we stress that these wavefunctions have nothing to do with quantum nature of the actual dark matter particles). We discuss different prescriptions to generate the initial wavefunctions from the initial conditions and demonstrate the validity of the technique on two simple test cases. This new simulation algorithm can in principle be used on an arbitrary distribution function, enabling the simulation of warm and hot dark matter structure formation scenarios.