Multiphase fluid flows are very common in engineering and science applications. Examples include air ow on water surface, metallurgical flow and blood flow in the body. In these flows, fluids are separated by a sharp interface and form different phases. The flow is characterized by the movement of this interface. Accurate modelling of the interface movement is a fundamental problem in the numerical simulation of these flows. Velocities for the movement are provided by the numerical solution of the Navier-Stokes (N-S) equations. These equations are discretized and converted into linear systems of equations. Research in the direction towards solving these systems efficiently has been the main focus of many researchers in the field of Computational Fluid Dynamics (CFD). A modified Volume of Fluid (VOF) method for modelling two phase flows is implemented using an analytic relation for its reconstruction step. The Finite Volume Method (FVM) is utilized, by incorporating a staggered grid, to discretize the two-dimensional (2-D) N-S equations. A preconditioned Krylov-Subspace iterative method, namely, the Bi-Conjugate Gradient Stabilized (Bi-CGSTAB) method is employed to solve the linear systems of equations. Solving the linear system usually consumes most of the simulation time for multiphase flow problems. Novel algorithms for the Incomplete LU Threshold (ILUT) preconditioner, forward and backward substitution and other matrix operations for penta-diagonal matrices are proposed here by adopting a diagonal sparse matrices format. The novel algorithm for ILUT reduces the computational complexity from O(n3 − n2) to O(n) in comparison to dense format. Further, it brings down the communication overhead, consequently facilitating parallelization. Parallel versions of these algorithms are developed using a new load balancing scheme. The MPI C++ communication library is utilized to develop the parallel version. The 2-D VOF code is applied to shape advection problems and results are found to be in good agreement with those available in literature. In the case of translation of a square box, it provides more accurate results than other VOF methods. The code for the VOF method and the parallel iterative solvers are integrated with 2-D N-S code in C++. The whole code is then implemented to simulate several two phase flow problems: dam breaking with and without an obstacle, rising of an air bubble and lid driven cavity flows. Speedup data from parallel programs implemented on these problems are generated.