Diffusion is the macroscopic manifestation of disordered molecular motion. Mathematically, diffusion equations are partial differential equations describing the fluid-like large-scale dynamics of parcels of molecules. Spatially inhomogeneous systems affect in a position-dependent way the average motion of molecules; thus, diffusion equations have to reflect somehow this fact within their structure. It is known since long that in this case an ambiguity arises: there are several ways of writing down diffusion equations containing space dependence within their parameters. These ways are all potentially valid but not equivalent, meaning that the different diffusion equations yield different solutions for the same data. The ambiguity can only be resolved at the microscopic level: a model for the stochastic dynamics of the individual molecules must be provided, and a well-defined diffusion equation then arises as the long-wavelength limit of this dynamics. In this work we introduce and employ the integro-differential Master Equation (ME) as a tool for describing the microscopic dynamics. We show that is possible to provide a parameterized version of the ME, in terms of a single numerical parameter (alpha), such that the different expressions for the diffusive fluxes are recovered at different values of alpha. This work aims to fill a gap in the literature, where the ME was shown to deliver just one diffusive limit. In the second part of the paper some numerical computer models are introduced, both to support analytical considerations, and to extend the scope of the ME to more sophisticated scenarios, beyond the simplest alpha-parameterization.