The accurate numerical solution of partial differential equations (PDEs) 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 PDEs 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. In contrast to previous works, this dynamical adaptation of the parameters is carried out using an explicit prescription avoiding iterative gradient descent. For the considered benchmark cases we observe excellent agreement with numerical solutions as well as analytical solutions for tasks that are challenging for traditional computational approaches.