Neville's algorithm is based on the Newton form of the interpolating polynomial and the recursion relation for the divided differences.
Let pi,j denote the polynomial of degree j − i which goes through the points (xk, yk) for k = i, i + 1, ..., j.
For instance, for n = 4, one can use the recurrence to fill the triangular tableau below from the left to the right.
This algorithm needs O(n2) floating point operations to interpolate a single point, and O(n3) floating point operations to interpolate a polynomial of degree n. The derivative of the polynomial can be obtained in the same manner, i.e: Lyness and Moler showed in 1966 that using undetermined coefficients for the polynomials in Neville's algorithm, one can compute the Maclaurin expansion of the final interpolating polynomial, which yields numerical approximations for the derivatives of the function at the origin.
They also show that their method can be applied directly to the solution of linear systems of the Vandermonde type.