The main objective of the current work is to enhance consistency and capabilities of Moving Particle Semi-implicit (MPS) method for simulating a wide range of free-surface flows and convection heat transfer. For this purpose, two novel high-order gradient and Laplacian operators are derived from the Taylor series expansion and are applied for the discretization of governing equations. Furthermore, the combination of the explicit Third-order TVD Runge-Kutta scheme and two-step projection algorithm is employed to approximate transient terms in the Navier-stokes and energy equations. To further improve the accuracy and performance of the method, a new kernel function is constructed by a combination of the Gaussian and cosine functions and then implemented for modeling the 1D Sod shock tube problem. Validation and verification of the proposed model are conducted through the simulations of several canonical test cases such as: dam break, rotation of a square patch of fluid, two-phase Rayleigh-Taylor instability, oscillating concentric circular drop and good agreement are achieved. The proposed model is then employed to simulate three-phase Rayleigh-Taylor instability and entropy generation due to natural convection heat transfer (Differentially Heated Cavity and Rayleigh-Bénard convection). The obtained results reveal that, the newly constructed kernel function provides more reliable results in comparison with two frequently used kernel functions namely; quartic spline and Wendland. Furthermore, it is found that, the enhanced MPS model is capable of handling multiphase flow problems with low and high density contrast.