We calculate the step scaling function, the lattice analog of the renormalization group β-function, for an SU(3) gauge theory with ten fundamental flavors. We present a detailed analysis including the study of systematic effects of our extensive data set generated with ten dynamical flavors using the Symanzik gauge action and three times stout smeared Möbius domain wall fermions. Using up to 32 4 volumes, we calculate renormalized couplings for different gradient flow schemes and determine the step-scaling β function for a scale change s = 2 on up to five different lattice volume pairs. In an accompanying paper we discuss that gradient flow can promote lattice dislocations to instantonlike objects, introducing nonperturbative lattice artifacts to the step scaling function. Motivated by the observation that Wilson flow sufficiently suppresses these artifacts, we choose Wilson flow with the Symanzik operator as our preferred analysis. We study systematic effects by calculating the step-scaling function based on alternative flows (Zeuthen or Symanzik), alternative operators (Wilson plaquette, clover), and also explore the effects of the perturbative tree-level improvement. Further we investigate the effects due to the finite value of Ls. c 11. At that strong coupling we also discovered previously unaccounted lattice artifacts. In an accompanying paper we discuss that gradient flow on coarse configurations can promote dislocations to instanton-like objects. This introduces a nonperturbative lattice artifact to the step scaling beta function which leads to incorrect continuum limit extrapolations [25]. We find that the perturbatively arXiv:2004.00754v1 [hep-lat] 2 Apr 2020 3 The good properties of this action has been first demonstrated for simulations in QCD [42,43] but is now also used in large scale simulations of a composite Higgs model [44][45][46]