function[b] = newton_int(N,x,f)

for k = 1:N
    P(k,1) = f(k);
end

for k = 2:N
    for j = 2:k
        P(k,j) = (P(k,j-1)-P(k-1,j-1))/(x(k)-x(k-j+1));
    end
end

for k = 1:N
    b(k) = P(k,k);
end