clc, clear, close all;
data = [1 3;2 8; 3 15];

x = -10:0.1:10;

y = 0;

for i = 1:size(data,1)
    p = ones(1,length(x));
    for j = 1:size(data,1)
        if i == j
            continue
        end
        p = p .* (x-data(j,1))/(data(i,1)-data(j,1));
    end
    y = y + p * data(i,2);
end

figure(1)
hold on;
plot(x,y);
plot(data(:,1), data(:,2), 'ko')