format long
x=[0.6 0.7 0.8 0.9];
c=[cos(0.6) cos(0.7) cos(0.8) cos(0.9)];
y=c-x;

n=size(x)(2)-1;
Q=zeros(n+1);
Q(:,1)=x';
z=0;
for i=1:n
  for j=1:i
    Q(i+1,j+1)=((z-y(i-j+1))*Q(i+1,j)-(z-y(i+1))*Q(i,j))/(y(i+1)-y(i-j+1));
  end
end
Q(n+1,n+1)