A Lagrangian particle multiphase model based on the Meshless Local Petrov-Galerkin method with Rankine source solution (MLPG_R) is proposed to simulate 2D flows of two immiscible fluids. The model is applicable to fluids with a wide range of density ratio from 1.01 to 1000, capable of dealing with violent flow situations (e.g. breaking waves) and maintaining the sharp discontinuity of properties of the fluids at the interface. In order to extend the MLPG_R method to model multiphase flows, an innovative phase coupling is proposed using an equation for pressure at the interface particle through two stages. In the first stage, the formulation is based on ensuring the continuity of the pressure and the ratio of the pressure derivative in the normal direction of the interface to the fluid density across the interface. Gravity current, natural sloshing of two layered liquids and air-water violent sloshing are successfully simulated and compared with either experimental data or analytical solution demonstrating a second order convergent rate. The second stage involves the extension of the method to account for interface tension and large viscosity. This is achieved by adding additional terms in the original pressure formulation considering jumps for both the pressure and the ratio of pressure derivatives in either the normal or tangential direction to fluid density at the interface. The method ensures both velocity continuity even with highly viscous fluids and interface stress balance in the presence of interface tension. Simulation of square-droplet deformation illustrates sharp pressure drop at the interface and relieved spurious currents that are known to be associated with predictions by other existing models. The capillary wave case also demonstrates the necessity of maintaining the jump in the ratio of pressure gradient to fluid density and the bubble rising case further validates the model as compared with the benchmark numerical results. Apart from being the first application of the MLPG_R method to multiphase flows the proposed model also contains two highly effective and robust techniques whose applicability is not restricted to the MLPG_R method. One is based on use of the absolute density gradient for identifying the interface and isolated particles which is essential to ensure that interface conditions are applied at the correct locations in violent flows. The effectiveness of the technique has been examined by a number of particle configurations, including those with different levels of randomness of particle distribution. The other is about solving the discretised pressure equation, by splitting the one set equations into two sets corresponding to two phases and solving them separately but coupled by the interface particles. This technique efficiently gives reasonable solutions not only for the cases with low density ratio but also for the ones with very high density ratio.