In infinite volume the gradient flow transformation can be interpreted as a continuous real-space Wilsonian renormalization group (RG) transformation. This approach allows one to determine the continuous RG $\beta$ function, an alternative to the finite-volume step-scaling function. Unlike step-scaling, where the lattice must provide the only scale, the continuous $\beta$ function can be used even in the confining regime where dimensional transmutation generates a physical scale $\Lambda_{\mathrm{QCD}}$. We investigate a pure gauge SU(3) Yang-Mills theory both in the deconfined and the confined phases and determine the continuous $\beta$ function in both. Our investigation is based on simulations done with the tree-level Symanzik gauge action on lattice volumes up to $32^4$ using both Wilson and Zeuthen gradient flow (GF) measurements. Our continuum GF $\beta$ function exhibits considerably slower running than the universal 2-loop perturbative prediction, and at strong couplings it runs even slower than the 1-loop prediction.