The Material Point Method (MPM), as proposed by Sulsky et al. (1994), has been developed to simulate large deformations and failure evolution involving different material phases in a single computational domain. A continuum body is divided into a finite number of subregions represented by Lagrangian material points, while the governing equations are formulated and solved with the Eulerian grid. Since this grid can be chosen arbitrarily, mesh tangling does not appear in the MPM. To design a simple but robust spatial discretization procedure, the MPM is coupled with the finite difference method (FDM) in the present study for simulating fully and partially saturated elasto-plastic soil responses based on the simplified three-phase method. Governing equations for the soil skeleton and the pore fluid are discretized by the MPM and FDM, respectively. Soil-water coupled analyses for fully saturated soils and seepage-deformation coupled analyses for unsaturated soils are performed, and the potential of the proposed method is demonstrated via numerical examples.