The accurate numerical solution of partial differential equations is a central task in numerical analysis allowing to model a wide range of natural phenomena by employing specialized solvers depending on the scenario of application. Here, we develop a variational approach for solving partial differential equations governing the evolution of high dimensional probability distributions. Our approach naturally works on the unbounded continuous domain and encodes the full probability density function through its variational parameters, which are adapted dynamically during the evolution to optimally reflect the dynamics of the density. For the considered benchmark cases we observe excellent agreement with numerical solutions as well as analytical solutions in regimes inaccessible to traditional computational approaches.