An integrated Eulerian-Lagrangian method has been proposed to simulate liquid jet in supersonic flow. The carrier fluid surrounding the liquid droplets is described by multicomponent Navier-Stokes equations based on stationary (Eulerian) Cartesian grid, and the Lagrangian-discrete droplet method is employed to represent the behavior of atomized droplets in this work. A total of three classic breakup models, Taylor analogy breakup model, Reitz wave model, and Kelvin-Helmholtz/Rayleigh-Taylor hybrid model, are discussed under supersonic conditions, and model predictions are compared to experimental data through multiple perspectives quantitatively. More accurate predictions of liquid penetration, as well as droplet size distribution, can be achieved for specific conditions with Kelvin-Helmholtz/Rayleigh-Taylor hybrid model comparing to other breakup models. Additionally, a statistical injection model has been introduced to depict droplet dispersion in the near-nozzle region. The probability density function is utilized to predict the size and position of the injected droplet, and the four commonly used probability density functions, uniform, chi-squared, Nukiyama-Tanasawa, and Rosin-Rammler, are analyzed. The influence of statistical representation of injection condition is discussed. The simulation results indicate that the random components in the velocity and droplet size of the injected droplets have a significant impact on the structure of the liquid jet in high-speed flow.