  We introduce a novel spatio-temporal discretization for nonlinear
Fokker-Planck equations on the multi-dimensional unit cube. This discretization
is based on two structural properties of these equations: the first is the
representation as a gradient flow of an entropy functional in the
$L^2$-Wasserstein metric, the second is the Lagrangian nature, meaning that
solutions can be written as the push forward transformation of the initial
density under suitable flow maps. The resulting numerical scheme is entropy
diminishing and mass conserving. Further, the scheme is weakly stable, which
allows us to prove convergence under certain regularity assumptions. Finally,
we present results from numerical experiments in space dimension $d=2$.
