  We perform the step-scaling investigation of the running coupling constant,
using the gradient-flow scheme, in SU(3) gauge theory with twelve massless
fermions in the fundamental representation. The Wilson plaquette gauge action
and massless unimproved staggered fermions are used in the simulations. Our
lattice data are prepared at high accuracy, such that the statistical error for
the renormalised coupling, g_GF, is at the subpercentage level. To investigate
the reliability of the continuum extrapolation, we employ two different lattice
discretisations to obtain g_GF. For our simulation setting, the corresponding
gauge-field averaging radius in the gradient flow has to be almost half of the
lattice size, in order to have this extrapolation under control. We can
determine the renormalisation group evolution of the coupling up to g^2_GF ~ 6,
before the onset of the bulk phase structure. In this infrared regime, the
running of the coupling is significantly slower than the two-loop perturbative
prediction, although we cannot draw definite conclusion regarding possible
infrared conformality of this theory. Furthermore, we comment on the issue
regarding the continuum extrapolation near an infrared fixed point. In addition
to adopting the fit ansatz a'la Symanzik for performing this task, we discuss a
possible alternative procedure inspired by properties derived from low-energy
scale invariance at strong coupling. Based on this procedure, we propose a
finite-size scaling method for the renormalised coupling as a means to search
for infrared fixed point. Using this method, it can be shown that the behaviour
of the theory around g^2_GF ~ 6 is still not governed by possible infrared
conformality.
