[en] We present a precise computation of the topological susceptibility χ_Y_M of SU(N) Yang-Mills theory in the large N limit. The computation is done on the lattice, using high-statistics Monte Carlo simulations with N=3,4,5,6 and three different lattice spacings. Two major improvements make it possible to go to finer lattice spacing and larger N compared to previous works. First, the topological charge is implemented through the gradient flow definition; and second, open boundary conditions in the time direction are employed in order to avoid the freezing of the topological charge. The results allow us to extrapolate the dimensionless quantity t_0"2χ_Y_M to the continuum and large N limits with confidence. The accuracy of the final result represents a new quality in the verification of large N scaling.