We develop here a rigorous homogenization approach applied to dense cubic arrays of magnetodielectric spheres, in order to determine their first-principle bulk constitutive parameters. We carefully model their electromagnetic response and assess the validity and limitations of commonly used approximate homogenization models for metamaterials. In particular, we show that, in regimes for which exotic wave propagation is obtained, near the array resonances and for especially dense arrays, more rigorous homogenization models are necessary, including magneto-electric coupling and weak spatial dispersion effects. We apply our findings to the design of negative-index metamaterial devices, and compare their scattering properties with those of our homogenized model. The comparison of near-field distributions shows remarkable agreement and confirms that a proper homogenization model can indeed correctly capture the physics and complex wave propagation in exotic metamaterial arrays.