We illustrate that a recently formulated recursive transfer matrix method can be used to reliably calculate the electromagnetic fields throughout three-dimensional systems of strongly scattering spheres, and/or coated spheres. The exceptional features of our technique are its particularly stable and reliable numeric implementations. In this work, we present new self-consistent formulae which permit the verification of the numerical stability at any stage of the calculations, and which ensure the satisfaction of the underlying multiple scattering equations for an arbitrary incident wave. Exact solutions to finite multiple scattering systems can provide important insights into the numerous approximate methods that are typically used to describe multiple scattering. Perhaps even more importantly, exact solutions are a good means of studying field fluctuations in multiple scattering systems.Invoking the addition theorem [1, 2], one can rather readily describe a multiple scattering system via coupled matrix equations acting on a spherical wave development of the incident field [3][4][5][6][7]. After a suitable truncation, the complete solution to this problem can then, in principle, be obtained by the inversion of an entire 'system' matrix [3,[5][6][7]. Analogous problems for two-dimensional systems of infinite cylinders have been treated with success in this manner [8,9]. The enormous number of unknowns and the sparse nature of the matrix, however, often prevents a direct matrix inversion solution for three-dimensional systems of spheres [5,6,[10][11][12].Several authors have proposed the calculation of cluster transfer matrices, T cl N , via recursive approaches which build the transfer matrix one particle at a time [11,13]. A major obstacle in these recursive approaches is that they suffer from numerical instabilities associated with partial-wave space truncations as recently discussed by Siqueira and Sarabandi [12], and by us [5]. Another deficiency of these new transfer-matrix approaches is that in practice, one often loses (or discards) local field information. In response to these difficulties, we have recently developed and applied a modified recursive technique which eliminates the aforementioned numerical difficulties [3][4][5].In this work, we illustrate that this new recursive transfermatrix technique can also be used to calculate the field distribution throughout multiple scattering systems (both near and far field). Since the numerical stability of the recursive transfer methods is still questioned in the multiple scattering community, we introduce here new self-consistence formulae which must be satisfied by our transfer matrices. We illustrate here that the satisfaction of these relations implies the satisfaction of the underlying multiple scattering equations for arbitrary incident wave configurations. To date, we have found that these relations are satisfied by our recursive matrix transfer matrices up to computer precision round off.