The group function theory by Tolpygo and McWeeny is a useful tool in treating quantum systems that can be represented as a set of localized electronic groups (e.g. atoms, molecules or bonds). It provides a general means of taking into account intra-correlation effects inside the groups without assuming that the interaction between the groups is weak. For non-orthogonal group functions the arrow diagram (AD) technique provides a convenient procedure for calculating matrix elements of arbitrary symmetrical operators which are needed, for example, for calculating the total energy of the system or its electron density. The total wavefunction of the system is represented as an antisymmetrized product of non-orthogonal electron group functions Φ_I of each group I in the system. However, application of the AD theory to extended (e.g. infinite) systems (such as biological molecules or crystals) is not straightforward, since the calculation of the mean value of an operator requires that each term of the diagram expansion be divided by the normalization integral S = Ψ|Ψ which is given by an AD expansion as well. In our previous work, we cast the mean value of a symmetrical operator in the form of an AD expansion which is a linear combination of linked (connected) ADs multiplied by numerical pre-factors. To obtain the pre-factors, a method based on power series expansion with respect to overlap was developed and tested for a simple 1D Hartree–Fock (HF) ring model. In the present paper this method is first tested on a 2D HF model, and we find that the power series expansion for the pre-factors converges extremely slowly to the exact solution. Instead, we suggest another, more powerful, method based on a continued fraction expansion of the pre-factors that approaches the exact solution much faster. The method is illustrated on the calculation of the electron density for the 2D HF model. It provides a powerful technique for treating extended systems consisting of a large number of strongly localized electronic groups.