This paper describes a parameter estimation algorithm applicable for the model structures in the form of multivariate polynomial systems. An example of a synchronous machine nonlinear model is used throughout the paper to explain the contribution. The fundamental in the proposed algorithm is the idea of reformulating the least squares estimation problem having the polynomial model into a numerical linear algebra problem. Recorded signals are represented using the Lagrange interpolation at the Chebyshev points which allows accurate computation of signal derivatives and numerical integration as both are required in the algorithm. The paper explores sensitivity of the algorithm to the power of recorded signals (i.e. excitation intensity) and discusses impact of roundoff error.