function Pu=prepare(X,Y)
n=length(X);
Pu=zeros(1,n);
P=zeros(n,n);

for i=2:n
    P(i,1)=Y(1);
end

for i=1:n-1
    for j=(i+2):n
        Pi=1;
        for k=1:i
            Pi=Pi*(X(j)-X(k))/(X(i+1)-X(k));
        end
        P(j,i+1)=P(j,i)+Pi*(Y(i+1)-P(i+1,i));
    end
end

for i=1:n-1
    Pu(i+1)=P(i+1,i);
end

return