f(x)=1/(x^2+1)
xs=[
 0.30567131760307564
 0.6509032415517444
 1.692835555726158
 2.665977427489991
 2.998121195753619
]
N=length(xs)
M=zeros(N,N)
for i=1:N 
    M[i,1]=f(xs[i])
end
for j=2:N
    for i=1:N-j+1
        M[i,j]=(M[i+1,j-1]-M[i,j-1])/(xs[i+j-1]-xs[i])
    end
end

function phi(n,x)
    r=1.0
    for k=1:n
        r*=(x-xs[k])
    end
    return r
end

function p(x,d)
    s=0
    for j=1:d+1
        s+=phi(j-1,x)*M[1,j]
    end
    return s
end