A new approach to molecular diffusion is developed using density functionals for fluxes and the Metropolis algorithm in the mass balance equation. This procedure results in a new equation for diffusion of interacting particles which has multiple solutions and gives density distributions for coexisting and metastable phases. It is shown that the diffusion of interacting molecules is driven by two variables: the density gradient of molecules and the density of molecule-vacancy pairs ("pseudoparticles").