A volume integral equation method (VIEM) is introduced for the solution of elastostatic problems in an unbounded isotropic elastic solid containing interacting multiple isotropic inclusions subject to uniform remote tension or in-plane shear. This method is applied to two-dimensional problems involving long parallel cylindrical inclusions. A detailed analysis of the stress field at the interface between the matrix, and the central inclusion is carried out for square and hexagonal packing of the inclusions. The effects of the number of isotropic inclusions and various fiber volume fractions on the stress field at the interface between the matrix and the central inclusion are also investigated in detail. The accuracy and efficiency of the method are examined in comparison with results obtained from analytical and finite element methods.