function p_xi = myNeville(x,y,xi)

aux=y;
n=length(x)-1;
for i=n:-1:1
    for j=1:i
        aux(j)=(xi-x(j))*aux(j+1)-(xi-x(j+n+1-i))*aux(j);
        aux(j)=aux(j)/(x(j+n+1-i)-x(j));
    end
end
p_xi=aux(1);