When multiple fluids in a single system interact, they can mix and separate according to their physical properties. We can describe a two-fluid system of this kind using the phase-field model. This model uses a scalar field to represent the volume concentration or concentration difference of one of the fluids. By varying this field we can describe phase mixtures, and the interface between phases is represented by the field smoothly transitioning between extreme values. The Cahn-Hilliard equation is a fourth-order partial differential equation that can be used to describe the change in this phase-field through separation. This method acts to minimize an energy potential function with minima where the mixture is fully separated. We couple the Cahn-Hilliard equation to the classic Navier-Stokes equations to include fluid flow. We develop two novel discretization methods for this system. The first method is a simplified time integration scheme based on splitting the fourth-order equation into a pair of second-order Helmholtz equations that can be solved as a pair of linear systems. We propose a method of coupling this discretization to an incompressible Navier-Stokes equation that maintains consistency between mass and momentum. We also propose a novel surface tension force discretization that limits surface forces to the interface region. This surface tension can be controlled with a parameter that allows for transitioning between a sharp interface formulation and a smooth continuum surface force. A drawback of this scheme is that it relies on a mass redistribution post-processing step to maintain the phase variable in the physical bounds. The second method is developed to intrinsically maintain the phase variable inside a physically consistent range instead of requiring a post-processing step. We formulate a piecewise energy potential function with energy barriers at both extremes. We combine this with a nonlinear second-order solve that is implicit in the energy potential. We show how to formulate this nonlinear equation as an optimization problem to be solved with a damped Newton's method. We compare both methods to previous work using a variety of experiments, and demonstrate that they are second-order accurate.