function [Pz,A]=interpol(x,y,z)
n=length(x);
A=zeros(n,n);
A(:,1)=y;
for j=2:n
    for i=j:n;
        A(i,j)=(A(i,j-1)-A(i-1,j-1))./(x(i)-x(i-j+1));
    end
end
Pz=A(n,n);
for k=n-1:-1:1
    Pz=A(k,k)+(z-x(k)).*Pz;
end
A=[x' ,A];
