Abstract: The numerical integration algorithm for the dual-yield surface constitutive model is crucial for accurately calculating the stress and deformation of soils. Based on a nonlinear elastoplastic model with a dual plastic mechanism, an implicit return mapping algorithm for the dual-yield surface model is established, and the corresponding consistent tangent modulus is derived. Considering the numerical singularities and convergence issues that can occur at the intersection of the dual yield surfaces when using the Newton algorithm for stress integration, a two-stage iterative algorithm is proposed in the plastic correction process. This involves first determining the initial iteration value using plastic increment theory, followed by iterative corrections using the Newton algorithm. Finally, a numerical solution program is developed using the UMAT interface provided by ABAQUS, and the model is analyzed and validated through triaxial test results. The results show that the numerical program effectively reflects the influence of different confining pressures on the stress-strain curves of sand, accurately describing the volumetric behavior of calcareous sand during dilatancy and compaction. Additionally, the new iterative scheme demonstrates superior computational efficiency compared to conventional iterative methods, effectively resolving numerical singularities and convergence issues, thereby highlighting the superiority, correctness, and practicality of the algorithm and the program.