We prove optimal convergence results of a stochastic particle method for computing the classical solution of a multivariate McKean-Vlasov equation, when the measure variable is in the drift, following the classical approach of [BT97, AKH02]. Our method builds upon adaptive nonparametric results in statistics that enable us to obtain a data-driven selection of the smoothing parameter in a kernel type estimator. In particular, we generalise the Bernstein inequality of [DMH21] for mean-field McKean-Vlasov models to interacting particles Euler schemes and obtain sharp deviation inequalities for the estimated classical solution. We complete our theoretical results with a systematic numerical study, and gather empirical evidence of the benefit of using high-order kernels and data-driven smoothing parameters.