We perform the step-scaling investigation of the running coupling constant,using the gradient-flow scheme, in SU(3) gauge theory with twelve masslessfermions in the fundamental representation. The Wilson plaquette gauge actionand massless unimproved staggered fermions are used in the simulations. Ourlattice data are prepared at high accuracy, such that the statistical error forthe renormalised coupling, g_GF, is at the subpercentage level. To investigatethe reliability of the continuum extrapolation, we employ two different latticediscretisations to obtain g_GF. For our simulation setting, the correspondinggauge-field averaging radius in the gradient flow has to be almost half of thelattice size, in order to have this extrapolation under control. We candetermine 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, therunning of the coupling is significantly slower than the two-loop perturbativeprediction, although we cannot draw definite conclusion regarding possibleinfrared conformality of this theory. Furthermore, we comment on the issueregarding the continuum extrapolation near an infrared fixed point. In additionto adopting the fit ansatz a'la Symanzik for performing this task, we discuss apossible alternative procedure inspired by properties derived from low-energyscale invariance at strong coupling. Based on this procedure, we propose afinite-size scaling method for the renormalised coupling as a means to searchfor infrared fixed point. Using this method, it can be shown that the behaviourof the theory around g^2_GF ~ 6 is still not governed by possible infraredconformality.