function lgip(order)

Nlg = order + 1;
Npts = 11;
t = Nlg/2 + ((0:Npts-1)/(Npts-1)-0.5)
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 * (t(k) - j)/((i)-(j));
            end;
        end;
        y(k) = y(k) + xk(i+1) * hlg;
    end;
end;
plot(y); grid;