Neuroimaging with magneto and electroencephalography (M/EEG) requires to compute the forward problem. It consists in predicting what is measured by MEG or EEG sensors due to a configuration of current generators within the head. When considering realistic head models, the equations derived from Maxwell equations can only be solved numerically. The Boundary Element Method is a standard approach to address this problem. However, different mathematical and computational variants exist within this class of BEM solvers. At least two BEM formulations exist for M/EEG forward modeling [Geselowitz 1967, Kybic et al. 2005]. Also implementation details vary between the different software packages (e.g. the precision of the numerical integrations). In order to investigate the influence of such differences, we have conducted a numerical experiment that evaluates the precision of different freely available BEM solvers. The BEM solvers are tested for EEG forward modeling with head models consisting of three nested shells. Multiple dipoles are located at different distances from the inner shell supposed to model the inner skull. To quantify the precision of the solvers we consider sphere models for which analytical solutions exist. Random sphere models are obtained by triangulating points randomly sampled over the spheres. We test and compare the accuracy of the BEM solvers on multiple random sphere models with different numbers of vertices on each interface. We also provide the computation times of the different solvers. The implementations tested are: OpenMEEG, SimBio [Zanow et al. 1995], Helsinki BEM [Stenroos et al. 2007], Dipoli [Oostendorp et al. 1989] and BEMCP [Phillips 2000] (Dipoli and BEMCP are available via Fieldtrip or SPM). Results show that OpenMEEG provides the best accuracy, followed by SimBio and dipoli which give comparable results. Finally, BEMCP and Helsinki BEM provide rather limited precisions.