This paper presents the recent progress achieved by the authors' group on the simulation of 2D elastic solids containing a large number of randomly distributed inclusions using BEM. A new scheme of repeated similar sub-domain BEM is proposed in this paper. The randomly distributed inclusions investigated include identical circular inclusions, elliptical inclusions with identical size and shape but different orientation, and also inclusions with different shape, size and different material properties. The number of inclusions is approximately 100 by conventional BEM, and it can be more than 1000 by the fast multipole BEM on only one PC. The interface between the matrix and inclusion can be at the ideal interface as well as at the interface with interphase layers. The numerical results show that for such investigations, the BEM is more suitable compared to FEM and other numerical methods. The simulated effective elastic modulus are presented and compared with different theoretical results. Future investigations should be extended to simulation of 3D models and simulation on brittle failure process.