The authors deal with the system of functional equations: \[ f_i (x) = \sum_{j=1}^n \sum_{k=1}^m a_{ijk} [x, f_j (s_{ijk} (x))] + g_i (x) , \] \(i= 1,2,...,n\) where \( x \in \Omega_i \), a compact or non-compact domain of \(\mathbb{R}^p \), \(g_i : \Omega_i \to \mathbb{R}\), \( S_{ijk} : \Omega_i \to \Omega_j \), \( a_{ijk} : \Omega_i \times \mathbb{R} \to \mathbb{R} \) are given continuous functions and \( f_i: \Omega_i \to \mathbb{R} \) are the unknown functions to be determined. The existence, uniqueness and stability of solutions are studied. Sufficient conditions to obtain quadratic convergence are given. Some particular cases are solved by means of Maclaurin expansions.