We present a numerical method for computing solutions of the incompressible Euler or Navier-Stokes equations when a principle feature of the flow is the presence of an interface between two fluids with different fluid properties such as density, viscosity, etc. The method is based on a second-order projection method for variable density flows using an ``approximate projection'' formulation. The boundary between the fluids is tracked with a second-order, volume-of-fluid interface tracking algorithm. We present results for a simple test problem with a known exact solution to demonstrate the accuracy of the interface tracking algorithm. We also present the results of computations of a Rayleigh-Taylor instability to demonstrate the behavior of the algorithm on a more difficult problem.