Due to the complex and dynamic nature of multi-phase interfaces, accurately capturing interface evolution remains one of the key challenges in multi-phase flow simulations, particularly in modeling bubble rising. In this study, a fully Lagrangian method is developed by using the Generalized Finite Difference (GFD) scheme, which we refer to as Finite Difference Particle Method (FDPM), and the Continuum Surface Force (CSF) model to simulate bubble dynamics. In this framework, the fluid is represented by particles, and all partial differential terms in the Navier–Stokes equations are discretized into symmetric linear systems using the GFD scheme. Notably, the interface curvature required by the CSF model is computed directly via the Laplacian of the color function, rather than through the divergence of the unit normal vector. The bubble relaxation cases with various density ratios (up to 1000) are tested, revealing that the pressure distributions inside and outside the bubbles generally agree with theoretical predictions. Finally, simulations of rising bubbles with a high density ratio (1:1000) and Reynolds number of 25 demonstrate that the time evolution of the bubble’s center of mass by FDPM is consistent with results obtained by the Smoothed Particle Hydrodynamics (SPH) and the Finite Element Method (FEM) approaches.