  We study the discrete beta function of SU(3) gauge theory with Nf=12 massless
fermions in the fundamental representation. Using an nHYP-smeared staggered
lattice action and an improved gradient flow running coupling $\tilde g_c^2(L)$
we determine the continuum-extrapolated discrete beta function up to $g_c^2
\approx 8.2$. We observe an IR fixed point at $g_{\star}^2 =
7.3\left(_{-2}^{+8}\right)$ in the $c = \sqrt{8t} / L = 0.25$ scheme, and
$g_{\star}^2 = 7.3\left(_{-3}^{+6}\right)$ with c=0.3, combining statistical
and systematic uncertainties in quadrature. The systematic effects we
investigate include the stability of the $(a / L) \to 0$ extrapolations, the
interpolation of $\tilde g_c^2(L)$ as a function of the bare coupling, the
improvement of the gradient flow running coupling, and the discretization of
the energy density. In an appendix we observe that the resulting systematic
errors increase dramatically upon combining smaller $c \lesssim 0.2$ with
smaller $L \leq 12$, leading to an IR fixed point at $g_{\star}^2 = 5.9(1.9)$
in the c=0.2 scheme, which resolves to $g_{\star}^2 =
6.9\left(_{-1}^{+6}\right)$ upon considering only $L \geq 16$. At the IR fixed
point we measure the leading irrelevant critical exponent to be
$\gamma_g^{\star} = 0.26(2)$, comparable to perturbative estimates.
