In this talk, we present numerical approximations of a phase-field model for two-phase ferrofluids, which consists of the Navier-Stokes equations, the Cahn-Hilliard equation, the magnetostatic equations, as well as the magnetic field equation. By combining the projection method for the Navier-Stokes equations and some subtle implicit-explicit treatments for coupled nonlinear terms, we construct a decoupled, linear, fully discrete finite element scheme to solve the highly nonlinear and coupled multi-physics system efficiently. The scheme is provably unconditionally energy stable and leads to a series of decoupled linear equations to solve at each time step. Through numerous numerical examples in simulating benchmark problems such as the Rosensweig instability and droplet deformation, we demonstrate the stability and accuracy of the numerical scheme.