We apply the Magnus expansion to the Zakharov-Shabat system, providing a basis for a systematic construction of high-order numerical schemes to solve the direct scattering problem of the integrable one-dimensional nonlinear Schrödinger equation. The presented numerical simulations of previously unreachable wave fields with up to 128 solitons employing second-, fourth-, and sixth-order schemes stresses the need for delicate numerics to identify the eigenvalues and especially phase coefficients. This approach lays the foundation for the study of large optical wave packets, providing fundamental information about their scattering data content and origin of various nonlinear effects.