  A model based on the classic non-interacting Ehrenfest urn model with
two-urns is generalized to $M$ urns with the introduction of interactions for
particles within the same urn. As the inter-particle interaction strength is
varied, phases of different levels of non-uniformity emerge and their
stabilities are calculated analytically. In particular, coexistence of locally
stable uniform and non-uniform phases connected by first-order transition
occurs. The phase transition threshold and energy barrier can be derived
exactly together with the phase diagram obtained analytically. These analytic
results are further confirmed by Monte Carlo simulations.
