  In the high energy limit of hadron collisions, the evolution of the gluon
density in the longitudinal momentum fraction can be deduced from the Balitsky
hierarchy of equations or, equivalently, from the nonlinear
Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation. The
solutions of the latter can be studied numerically by using its reformulation
in terms of a Langevin equation. In this paper, we present a comprehensive
study of systematic effects associated with the numerical framework, in
particular the ones related to the inclusion of the running coupling. We
consider three proposed ways in which the running of the coupling constant can
be included: "square root" and "noise" prescriptions and the recent proposal by
Hatta and Iancu. We implement them both in position and momentum spaces and we
investigate and quantify the differences in the resulting evolved gluon
distributions. We find that the systematic differences associated with the
implementation technicalities can be of a similar magnitude as differences in
running coupling prescriptions in some cases, or much smaller in other cases.
