We present a precise computation of the topological susceptibility $\chi_{_\mathrm{YM}}$ 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\chi_{_\mathrm{YM}}$ 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.