we develop and analyze a numerical method proposed for solving high-dimensional fokker-planck equations by leveraging the generative models from deep learning. our starting point is a formulation of the fokker-planck equation as a system of ordinary differential equations (odes) on finite-dimensional parameter space with the parameters inherited from generative models such as normalizing flows. we call such odes "neural parametric fokker-planck equations". the fact that the fokker-planck equation can be viewed as the 2-wasserstein gradient flow of the relative entropy (also known as kl divergence) allows us to derive the ode as the 2-wasserstein gradient flow of the relative entropy constrained on the manifold of probability densities generated by neural networks. for numerical computation, we design a bi-level minimization scheme for the time discretization of the proposed ode. such an algorithm is sampling-based, which can readily handle computations in higher-dimensional space. moreover, we establish bounds for the asymptotic convergence analysis as well as the error analysis for both the continuous and discrete schemes of the neural parametric fokker-planck equation. several numerical examples are provided to illustrate the performance of the proposed algorithms and analysis.