Abstract: This contribution is the numerically oriented companion article of the work [E. Canc\`es, V. Ehrlacher, F. Legoll, B. Stamm and S. Xiang, arxiv preprint 1807.05131]. We focus here on the numerical resolution of the embedded corrector problem introduced in [E. Canc\`es, V. Ehrlacher, F. Legoll and B. Stamm, CRAS 2015; E. Canc\`es, V. Ehrlacher, F. Legoll, B. Stamm and S. Xiang, arxiv preprint 1807.05131] in the context of homogenization of diffusion equations. Our approach consists in considering a corrector-type problem, posed on the whole space, but with a diffusion matrix which is constant outside some bounded domain. In [E. Canc\`es, V. Ehrlacher, F. Legoll, B. Stamm and S. Xiang, arxiv preprint 1807.05131], we have shown how to define three approximate homogenized diffusion coefficients on the basis of the embedded corrector problems. We have also proved that these approximations all converge to the exact homogenized coefficients when the size of the bounded domain increases. We show here that, under the assumption that the diffusion matrix is piecewise constant, the corrector problem to solve can be recast as an integral equation. In case of spherical inclusions with isotropic materials, we explain how to efficiently discretize this integral equation using spherical harmonics, and how to use the fast multipole method (FMM) to compute the resulting matrix-vector products at a cost which scales only linearly with respect to the number of inclusions. Numerical tests illustrate the performance of our approach in various settings.