Numerical modeling of rock-breaking by a waterjet is a complex problem involving various nonlinear and coupled phenomena such as high velocity, large strains and fluid–structure interaction. This paper introduces the coupled Eulerian–Lagrangian (CEL) method, in which the water is described by the Eulerian formulation and the rock by the Lagrangian formulation. Furthermore, it adopts the modified Drucker–Prager (D-P) constitutive model to consider the effects of strain hardening, damage evolvement of the rock material, and rational failure criterion. This paper investigates three jet velocities to analyze the formation of the breaking hole, the distribution of stress, the depth of the broken hole, and the shape of the water flow, which are generally consistent with the experimental and numerical data in the literature. This result indicates that the application of the CEL method provides the possibility of modeling geotechnical problems of rock-breaking by a high-pressure waterjet.