To formulate estimation problem for a observed nonlinear system, differential equations are required as a system model. We recast those equations as integral equations and apply integral operator, which is approximated by Chebyshev interpolating polynomial. Resulting equations of the nonlinear estimation problem are integral equations discretized on the grid of Chebyshev points. A sequence of linear least squares problems are solved iteratively. Grid resolution can be determined automatically to maximize computation accuracy. Numerical efficiency is achieved by applying iterative method that requires only matrix-vector multiplications, and via implementation of Discrete Cosine Transform when solving indefinite integrals. The estimation algorithm works with matrices having bounded low condition number compared to large and unbounded condition number for the formulation with differential operator. This achievement has important practical value when applying the algorithm with high-order models when the differential operator formulation is typically ill-conditioned. Application to a Duffing system having chaotic response, has been used to illustrate advantages of the proposed estimation algorithm based on Chebyshev integral operator.