In this study, a higher order compact (HOC) finite difference scheme is proposed for the simulation of multiphase flows with surface tension effects by coupling the pure streamfunction formulation of Navier–Stokes (N–S) equations with the level set method. Unlike most of the existing methods for multiphase flows that solve the incompressible Navier–Stokes equations in the velocity–pressure form or the stream-function vorticity form, we recast the governing equations into the biharmonic form of the N–S equations which carries stream-function $\psi$ as the only variable to be computed. The coupling of this formulation with level set (LS) enables the simulation of two phase flows by the core computation of only the single variable $\psi$ free from discontinuities in the whole physical domain. In the process, we also propose a novel hybrid HOC algorithm for post-processing the pressure on either side of the interface including the jump across it. The robustness of this method is established through a series of benchmark problems such as dispersion of capillary waves, static bubble and oscillating bubble at equilibrium in a fluid at rest and deforming circular drops in shear flows. We compare our numerical results with available benchmark results and excellent match is obtained in all the cases along with accurate representation of overall trends. Our approach is proved to be very efficient in terms of its easy implementation, relatively low computational time and high order of accuracy.