Newton's interpolation

function p = newton_interpolation(x,y)

    cols = length(x);
    p = zeros(1, cols);
    D = zeros(cols, cols);
    D(1,:) = y';

    for j = 2:cols
        for i = 1:n-j+1
            D(i, j) = (D(i+1,j-1) - D(i,j-1))/(x(i+j-1) - x(i));
        end
    end

    for i=1:n
        pom = 1;
        for j=2:i
            pom = conv(pom,[1 -x(j-1)]);
        end
        p = p + [zeros(1,n-i) pom * D(1,i)];
    end

end