A volume integral equation method is introduced for the solution of elastostatic problems in an unbounded isotropic elastic solid containing interacting multiple anisotropic elliptical inclusions subject to uniform remote tension or remote in-plane shear. This method is applied to two-dimensional problems involving long parallel elliptical cylindrical inclusions. A detailed analysis of the stress field at the interface between the matrix and the central elliptical inclusion is carried out for square and hexagonal packing of anisotropic inclusions. The effects of the number of anisotropic inclusions and various inclusion volume fractions on the stress field at the interface between the isotropic matrix and the central elliptical cylindrical inclusion are investigated in detail. The stress field at the interface between the isotropic matrix and the central elliptical inclusion is also compared with that between the isotropic matrix and the central circular inclusion.