We study the system of functional equations f_i(x) = \sum^n_{j=1} \sum^m_{k=1} a_{ijk} [x, f_i (S_{ijk} (x))] + g_i (x) \ \ (l ≤ i ≤n) for x \in \Omega_i where \Omega_i are compact or non-compact domains of \mathbb R^p, g_i : \Omega_i \to 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 unknown functions. The paper consists of two mains parts. In the first part we give some results on existence, uniqueness and stability of the solutions of such systems and study sufficient conditions to obtain quadratic convergence. In the second part we obtain the Maclaurin expansion and approximation of solution in the case that a_{ijk} are linear and S_{ijk} are affine functions.