A Multiscale Finite Volume (MSFV) method is provided to efficiently solve large heterogeneous problems; it is usually employed for pressure equations and delivers conservative flux fields to be used in transport problems. It relies on the hypothesis that the fine-scale problem can be described by a set of local solutions coupled by a conservative coarse-scale problem. In numerically challenging cases, a more accurate localization approximation is used to obtain a good approximation of the fine-scale solution. According to an embodiment, a method is provided to iteratively improve the boundary conditions of the local problems, and is responsive to the data structure of the underlying MSFV method and employs a Krylov-subspace projection method to obtain an unconditionally stable scheme and accelerate convergence. In one embodiment the MSFV operator is used. Alternately, the MSFV operator is combined with an operator derived from the problem solved to construct the conservative flux field.