We calculate Eshelby tensor for inclusions of non-ellipsoidal shape. We focus on the superspherical shape described by equation X-2P + y(2p) + Z(2p) <= 1. It is convex when p > 0.5 and concave when p < 0.5. We propose a numerical approach to perform: integration on the surface of the superspherical inclusion necessary to compute the average Eshelby tensor. Validation of the method is done by comparison of the results with analytical solutions for a spherical inclusion (p = 1) and with numerical results of Onaka (2001) (p > 1).