f=@(x) cos(x);
x=[0.6 0.7 0.8 0.9];
y=[f(x(1)) f(x(2)) f(x(3)) f(x(4))]-x;
z=0;
q=zeros(4,4);
q(1,1)=x(1);
q(2,1)=x(2);
q(3,1)=x(3);
q(4,1)=x(4);
for i=1:1:3
  for j=1: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
    printf('%.10d',q(4,4));