This paper discusses approaches to the numerical integration of the coupled nonlinear Schrödinger equations system, different from the generally accepted approach based on the method of splitting according to physical processes. A combined explicit/implicit finite-difference integration scheme based on the implicit Crank–Nicolson finite-difference scheme is proposed and substantiated. It allows the integration of a nonlinear system of equations with a choice of nonlinear terms from the previous integration step. The main advantages of the proposed method are: its absolute stability through the use of an implicit finite-difference integration scheme and an integrated mechanism for refining the numerical solution at each step; integration with automatic step selection; performance gains (or resolutions) up to three or more orders of magnitude due to the fact that there is no need to produce direct and inverse Fourier transforms at each integration step, as is required in the method of splitting according to physical processes. An additional advantage of the proposed method is the ability to calculate the interaction with an arbitrary number of propagation modes in the fiber.