We consider accurate and iteratively efficient solutions of electromagnetic problems involving homogenized near-zero-index (NZI) bodies using surface-integral-equation formulations in the frequency domain. NZI structures can be practically useful in a plethora of optical applications, as they possess near-zero permittivity and/or permeability values that cannot be found in nature. Hence, numerical simulations are of the utmost importance for rigorous design and analysis of NZI structures. Unfortunately, small values of electromagnetic parameters bring computational challenges in numerical solutions of homogeneous models. Conventional formulations available in the literature encounter stability issues that make them inaccurate and/or inefficient as permittivity and/or permeability approach zero. We propose a novel formulation that involves a well-balanced combination of operators and that can provide both accurate and efficient solutions for all NZI cases. Numerical results are presented to demonstrate the superior properties of the developed formulation in comparison to the conventional ones.