The aim of this work is to develop a general numerical technique for obtaining the equilibrium shape of a precipitate in a matrix, taking into consideration the anisotropy in both interfacial and elastic bulk energies. Such shape optimization problems in previous literature have altogether ignored the anisotropy of the interfacial energy. The computational technique used in this work is based on the combination of an Extended finite element method (XFEM) and a level set method (LSM). The XFEM technique calculates the stress field in the system, with a moving interface, and the level set method updates the position of the interface on the basis of the velocity field obtained from the interfacial and elastic components. An isolated precipitate as well as the interaction between multiple precipitates are studied. After verifying our framework, the problems are revisited taking anisotropy in the interfacial energy into consideration. The issues of stability, bifurcation, symmetry loss, and interaction are discussed in detail.