Ellipsoidal inclusions can be regarded as a general model of defects in structures because they cover a lot of particular cases, such as line, circular and spherical defects. This paper deals with three-dimensional stress analysis for ellipsoidal inclusions in a bimaterial body under tension. The problem is formulated as a system of singular equations with Cauchy-type or logarithmic-type singularities, where unknowns are densities of body forces distributed in the r-and z-directions in bimaterial bodies having the same elastic constants of those of the given problem. In order to satisfy the boundary conditions along the ellipsoidal boundary, four types of fundamental density functions proposed in the previous paper are applied. Then the body force densities are approximated as a linear combination of fundamental density functions and polynomials. The present method is found to yield rapidly converging numerical results for stress distribution along the boundaries of both the matrix and inclusion even when the inclusion is very close to the bimaterial interface. Then, the effect of bimaterial surface on the stress concentration factor is discussed with varying the distance from bimaterial interface, shape ratio, and elastic modulus ratio.