function lgip(order)

Nlg = order + 1;
Npts = 100;
t = (0:Npts-1)/Npts;
xk = [0 0.5 0];
for k=1:Npts,
    y(k) = 0;
    for i=0:Nlg-1,
        hlg = 1;
        for j=0:Nlg-1,
            if (i ~= j)
                hlg = hlg * (Nlg/2 - 0.5 + t(k) - j)/(i-j);
            end;
        end;
        y(k) = y(k) + xk(i+1) * hlg;
    end;
end;
plot(t, y); grid;