This paper discusses approaches to the numerical integration of the coupled nonlinear Schrödinger equations system in case of few-mode wave propagation. The wave propagation assumes the propagation of up to nine modes of light in an optical fiber. In this case, the light propagation is described by the non-linear coupled Schrödinger equation system, where propagation of each mode is described by own Schrödinger equation with other modes interactions. In this case, the non-linear coupled Schrödinger equation system solving becomes increasingly complex, because each mode affects the propagation of other modes. The suggested solution is based on the direct numerical integration approach, which is based on a finite-difference integration scheme. The well-known explicit finite-difference integration scheme approach fails, due to the non-stability of the computing scheme. Due to this fact, the combined explicit/implicit finite-difference integration scheme, based on the implicit Crank–Nicolson finite-difference scheme, is used. It allows ensuring the stability of the computing scheme. Moreover, this approach allows separating the whole equation system on the independent equation system for each wave mode at each integration step. Additionally, the algorithm of numerical solution refining at each step and the integration method with automatic integration step selection are used. The suggested approach has performance gains (or resolutions) up to three or more orders of magnitude in comparison with the split-step Fourier method due to the fact that there is no need to produce direct and inverse Fourier transforms at each integration step. The main advantage of the proposed method is the ability to calculate the propagation of an arbitrary number of modes in the fiber.