The flow of two immiscible fluids in porous media is ubiquitous particularly in petroleum exploration and extraction. The displacement of one fluid by another immiscible with it represents a very important aspect in what is called enhanced oil recovery. Another example is related to the long-term sequestration of carbon dioxide, CO2 , in deep geologic formations. In this technique, supercritical CO2 is introduced into deep saline aquifer where it displaces the hosting fluid. Furthermore, very important classes of contaminants that are very slightly soluble in water and represent a huge concern if they get introduced to groundwater could basically be assumed immiscible. These are called light non-aqueous phase liquids (LNAPL) and dense non-aqueous phase liquids (DNAPL). All these applications necessitate that efficient algorithms be developed for the numerical solution of these problems. In this work we introduce the use of shifting matrices to numerically solving the problem of two-phase immiscible flows in the subsurface. We implement the cell-center finite difference method which discretizes the governing set of partial differential equations in conservative manner. Unlike traditional solution methodologies, which are based on performing the discretization on a generic cell and solve for all the cells within a loop, in this technique, the cell center information for all the cells are obtained all at once without loops using matrix oriented operations. This technique is significantly faster than the traditional looping algorithms, particularly for larger systems when coding using languages that require repeating interpretation each time a loop is called like Mat Lab, Python and the like. We apply this technique to the transport of LNAPL and DNAPL into a rectangular domain.