  Nonperturbative determinations of the renormalization group (RG) $\beta$
function are crucial to understand properties of gauge-fermion systems at
strong coupling and connect lattice simulations and the perturbative
ultraviolet regime. Choosing well-understood, QCD-like systems with SU(3) gauge
group and either six or four fundamental flavors, we investigate their
step-scaling $\beta$ function. In both cases we push the simulations to the
boundary of chiral symmetry breaking and study the regime $g^2_{GF} \lesssim
8.2$ with six, and $g^2_{GF} \lesssim 6.6$ with four flavors. We carefully
consider the lattice discretization errors by comparing three different
gradient flows (GF), and for each flow three operators to estimate the
renormalized finite volume coupling. We also consider the tree level
improvement of the coupling. Noteworthy outcome is that nonperturbatively
determined $\beta$ functions run much slower than perturbatively predicted.
