  We investigate the precision of the numerical implementation of the
functional renormalization group based on extracting the eigenvalues from the
linearized RG transformation. For this purpose, we implement the LPA and
$O(\partial^2)$ orders of the derivative expansion for the three-dimensional
$O(N)$ models with $N~\in~\{1,2,3\}$. We identify several categories of
numerical error and devise simple tests to track their magnitude as functions
of numerical parameters. Our numerical schemes converge properly and are
characterized by errors of several orders of magnitude smaller than the error
bars of the derivative expansion for these models. We highlight situations in
which our methods cease to converge, most often due to rounding errors. In
particular, we observe an impaired convergence of the discretization scheme
when the $\tilde \rho$ grid is cut off at the value $\tilde \rho_{\text{Max}}$
smaller than $3.5$ times the local potential minimum. The program performing
the numerical calculations for this study is shared as an open-source library
accessible for review and reuse.
