The flow of two or more immiscible fluids in porous media is widespread, particularly in the oil industry. This includes secondary and tertiary oil recovery and carbon dioxide (CO2) sequestration. Accurate predictions of the development of these processes are important in estimating the benefits and consequences of the use of certain technologies. However, this accurate prediction depends—to a large extent—on two things. The first is related to our ability to correctly characterize the reservoir with all its complexities; the second depends on our ability to develop robust techniques that solve the governing equations efficiently and accurately. In this work, we introduce a new robust and efficient numerical technique for solving the conservation laws that govern the movement of two immiscible fluids in the subsurface. As an example, this work is applied to the problem of CO2 sequestration in deep saline aquifers; however, it can also be extended to incorporate more scenarios. The traditional solution algorithms to this problem are modeled after discretizing the governing laws on a generic cell and then proceed to the other cells within loops. Therefore, it is expected that calling and iterating these loops multiple times can take a significant amount of computer time. Furthermore, if this process is performed with programming languages that require repeated interpretation each time a loop is called, such as Matlab, Python, and others, much longer time is expected, particularly for larger systems. In this new algorithm, the solution is performed for all the nodes at once and not within loops. The solution methodology involves manipulating all the variables as column vectors. By use of shifting matrices, these vectors are shifted in such a way that subtracting relevant vectors produces the corresponding difference algorithm. It has been found that this technique significantly reduces the amount of central-processing-unit (CPU) time compared with a traditional technique implemented within the framework of Matlab.