This paper deals with numerical solutions of singular integral equations in interaction problems of elliptical inclusions under general loading conditions. The stress and displacement fields due to a point force in infinite plates are used as fundamental solutions. Then, the problems are formulated as a system of singular integral equations with Cauchy-type or logarithmic-type singularities, where the unknowns are the body force densities distributed in infinite plates having the same elastic constants as those of the matrix and inclusions. To determine the unknown body force densities to satisfy the boundary conditions, four auxiliary unknown functions are derived from each body force density. It is found that determining these four auxiliary functions in the range 0≦φk≦π/2 is equivalent to determining an original unknown density in the range 0≦φk≦2π. Then, these auxiliary unknowns are approximated by using fundamental densities and polynomials. Initially, the convergence of the results such as unknown densities and interface stresses are confirmed with increasing collocation points. Also, the accuracy is verified by examining the boundary conditions and relations between interface stresses and displacements. Randomly or regularly distributed elliptical inclusions can be treated by combining both solutions for remote tension and shear shown in this study.