  We develop a model describing the behavior of two-phase ferrofluid flows
using phase field-techniques and present an energy-stable numerical scheme for
it. For a simplified, yet physically realistic, version of this model and the
corresponding numerical scheme we prove, in addition to stability, convergence
and as by-product existence of solutions. With a series of numerical
experiments we illustrate the potential of these simple models and their
ability to capture basic phenomenological features of ferrofluids such as the
Rosensweig instability.
