We propose an analytical solution of the multi-density Ornstein-Zernike equation supplemented by the associative Percus-Yevick closure relations specifically designed to describe the equilibrium properties of the novel class of patchy colloidal particles represented by the inverse patchy colloids with arbitrary number of patches. Using Baxter's factorization method, we reduce solution of the problem to the solution of one nonlinear algebraic equation for the fraction of the particles with one non-bonded patch. We present closed-form expressions for the structure (structure factor) and thermodynamic (internal energy) properties of the system in terms of this fraction (and parameters of the model). We perform computer simulation studies and compare theoretical and computer simulation predictions for the pair distribution function, internal energy, and number of single and double bonds formed in the system, for two versions of the model, each with two and three patches. We consider the models with formation of the double bonds blocked by the patch-patch repulsion and the models without patch-patch repulsion. In general very good agreement between theoretical and computer simulation results is observed.