A model has been developed to determine the effective ac magnetic response of magnetic systems, taking into account the demagnetization effects arising from the sample geometry which determine the out-of-phase components of the applied fundamental frequency and higher harmonic components. Indeed, demagnetization fields and their intermodulation can significantly affect the ac magnetic response. This approach provides a system of self-consistent linear equations relating the magnetic response to the external magnetic field by means of nonlinear magnetic susceptibility. The model is extended to the magnetic response of granular systems in terms of the contributions of the individual grains and of the whole sample in the presence of demagnetization effects of the whole sample and of the grains on a macroscopic scale. In particular, our model is applied to a granular superconducting system. The comparison between the performed numerical simulations and the experimental data shows that the demagnetization fields of the single grains and of the whole sample, and their intermodulation, are relevant if magnetic measurements are used to extract detailed information about the analyzed material.