A volume integral equation method (VIEM) is used to calculate the plane elastostatic field in an isotropic elastic half-plane containing multiple isotropic or anisotropic inclusions subject to remote loading. A detailed analysis of stress field at the interface between the matrix and the central inclusion in the first column of square packing is carried out for different values of the distance between the center of the central inclusion in the first column of square packing of inclusions and the traction-free surface boundary in an isotropic elastic half-plane containing multiple isotropic or anisotropic inclusions. The method is shown to be very accurate and effective for investigating the local stresses in an isotropic elastic half-plane containing multiple isotropic or anisotropic inclusions.