In this paper, a robust interfacial lattice Boltzmann flux solver (ILBFS) is presented to simulate multiphase flows, which is able to deal with large density contrasts. Different from the conventional finite-volume multiphase flow solvers, which calculate the fluxes at each cell interface based on macroscopic variables, the present method employs the Chapman-Enskog analysis to re-construct the convective and diffusive terms of the macroscopic governing equation, and uses the local equilibrium solutions of the lattice Boltzmann equation (LBE) for calculating the fluxes through the cell interface. As the phase information of multiphase flows is naturally brought by the LBE from both sides of the interface, the ILBFS simulating interfacial phenomena is more physical and precise. In addition, different from the conventional LB model for simulating the multiphase flows, the proposed ILBFS only applies LBE locally at each cell interface. Macroscopic variables at cell centers are given by directly solving the governing partial differential equations. Moreover, the physical boundary conditions in the present solver can be directly implemented by using the macroscopic variables, and are not necessary to transform to the LB distribution functions. Hence, the current solver requires less memory resources for data storage as only the equilibrium distribution function related to the conservative variables is involved. To validate the present solver, several benchmark cases, including the Laplace law, two stationary bubbles without collision, Rayleigh-Taylor instability, bubble rising under buoyancy and droplet splashing on a thin film at density radio of 1000, are studied. The obtained numerical results are in good agreement with the analytical solutions or the published results in the literature. (C) 2019 Elsevier Ltd. All rights reserved.