f=@(x)cos(x);
x=[0.6,0.7,0.8,0.9];
y(:,1)=f(x)-x;
y=y';
p=0;
k=0;
l=1;
for i=1:4
  for j=1:4
    if i!=j
      l=l*(k-y(j))/(y(i)-y(j));
    end
  end
  p=p+l*x(i);
  l=1;
end