A model based on the classic noninteracting 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 nonuniformity emerge and their stabilities are calculated analytically. In particular, coexistence of locally stable uniform and nonuniform 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.