Função para Cálculo pelo Método Polinômio Interpolador de Newton para Gretl

function coef=polinomio_interp_newton(x,y)

    // coef=polinomio_interp_newton(x,y)

    // Calcula os coeficientes do polinômio interpolador de Newton

    m = max(size(x));

    n = m-1;

    c = zeros(m,m);

    c(:,1) = y';

    for j=1:n

        for i=0:n-j

            c(i+1,j+1) = (c(i+2,j)-c(i+1,j))/(x(i+j+1)-x(i+1));

        end

    end

    coef = c(1,1:m);

endfunction