pkg load symbolic

xx = [0, 1, 3, 4];
yy = [1, 0, 2, 1];

s = 0;

syms x

for i = 1:length(xx)
    p = 1;
    for k = [1:i-1, i+1:length(xx)]
        p = p * (x - xx(k))/(xx(i) - xx(k));
    endfor

    s = s + p*yy(i);

endfor

s