The commonly used Eulerian or continuum model for incompressible multiphase flow is known to be unstable to perturbations for all wavenumbers, even if viscosity terms are used in the momentum equations. In the present work the model is stabilized by adding explicit artificial diffusion to the mass equations. The artificial diffusion terms lead to improved stability properties: uniform flow becomes linearly stable for large wavenumbers, and above an analytically derived threshold for the artificial diffusivity, stability for all wavenumbers is achieved. The artificial diffusivity reappears in the momentum equations, in such a way that fundamental properties of the standard equations remain valid: Galilean invariance is maintained, total mass and momentum are conserved, decay of total kinetic energy is ensured in the absence of external forces, and a flow initially at rest at hydrostatic pressure remains unchanged, even if the spatial distribution of volume fractions is nonuniform. A staggered finite volume pressure correction method using central differencing (leading to energy conserving discretization of convective and pressure terms) is presented. Application of the method to one-dimensional two-phase flow of falling particles particles confirms that the equations are stable with and unstable without artificial diffusion in the volume fraction equation.