procedura Lagrange()

> Lagrange := proc(k,body)
        global L;
        local a,i,j,n;

        a:=array(1..k):
        for i from 1 by 1 to k do
            a[i]:= product( (x-body[j][1])/(body[i][1]-body[j][1]), j=1..i-1 )*
            product( (x-body[j][1])/(body[i][1]-body[j][1]), j=i+1..k )
        od;
        L(x):=sum(body[n][2]*a[n], n=1..k);
        L(x) := expand(L(x));
    end: