This paper presents an iterative method for numerical solution of a system of n nonlinear functional equations in n unknowns, by repeated linear interpolation. It is a variation of Newton’s method, in which the partial derivatives are replaced by the multidimensional analogues of difference quotients. The only computations required are evaluation of the given functions and simple matrix operations. The method is a generalization of an idea of C. F. Gauss [6], and can also be regarded as generalizing the secant method, or regula falsi, as used for solving functions of one variable [4]. It produces superlinear convergence (of order $\frac{1}{2}(1 + \sqrt 5 )$, or approximately 1.62) to a simple zero, provided the given functions are twice continuously differentiable in a neighborhood of the zero, and the initial approximations are close enough. A convergence analysis for this method is given, and some numerical results are presented.