A combined numerical and analytic approach is presented to obtain a general matrix form of constitutive relations for a bi-anisotropic material which consists of arbitrary inclusions in a host medium. This approach is based on the quasi-static Lorentz-Lorenz theory which relates the constitutive parameters of the material with the electric and magnetic dipole moments of the inclusions which can be calculated in a numerical or analytic way. Analysis of "meta-materials" is conducted to validate the proposed method.