format long

d = 3;
n = d;
z = 0;

x = [0.6 0.7 0.8 0.9];
fx = cos(x);
fx = fx-x;

for k = 1:n+1
  L(k) = 1;
  for i = 1:n+1
    if(i~=k)
      L(k) = L(k)*(z-fx(i))/(fx(k) - fx(i));
    end
  end
end

ans = 0;

for a = 1:n+1
  ans = ans + x(a)*L(a);
end

ans