Fundamental solutions of elastostatics for infinite anisotropic media are obtained by numerical integration, and also by the finite element method for bounded sphere. These solutions are presented in the form of mean-squared approximation by the series of spherical harmonics Ylm up to the order l = 10 (look for Online Support Data, pdf) as exemplified by the materials of different elastic symmetry: isotropic (concrete), cubic (Si), hexagonal (AlN), orthorhombic (MgSiO3 in perovskite phase), tetragonal with 6 or 7 independent elastic constants (ZrSiO4 and CaWO4 respectively), trigonal with 6 or 7 constants (Al2O3 and dolomite), monoclinic (gypsum), triclinic (Al2SiO5). Using obtained Green's functions for each crystal, the energy of elastic interaction of a pair of point defects has been plotted as a function of angles of their mutual orientation.