Příklady

> A:=array(1..k,1..k):
    for i from 1 by 1 to k do
        A[i,1]:=1:
    end do:
    for i from 1 by 1 to k do
        for j from 2 by 1 to k do
            A[i,j]:=0:
        end do:
    end do:
    f(y):=1:
    for j from 2 by 1 to k do
        f(y):=f(y)*(y-body[j-1][1]):
        for i from j by 1 to k do
            A[i,j]:=eval(f(y),y=body[i][1]):
        end do:
    end do:
    eval(A);

> B:=array(1..k):
    for i from 1 by 1 to k do
        B[i]:=body[i][2]:
    end do:
    eval(B);

... a vyřešíme ji.

> a:=linsolve(A,B);

Nelezené řešení, tedy koeficienty a₁ ... a₄ už můžeme dosadit do vzorce ...

> N(x):=sum(a[n]*product((x-body[m][1]),m=1..n-1), n=1..k);

... a zjednodušit.

> N(x):=expand(N(x));

Graf opravdu prochází všemi čtyřmi zadanými uzlovými body.

> plots[display]({plot(body,x=body[1][1]..body[k][1],style=point,color=blue),
                            plot(N(x),x=body[1][1]..body[k][1],style=line,color=red)
                           });

2. příklad:

Vytvořme si tabulku funkce f(x)  = sin(x) na intervalu <-Pi ; Pi> tak, že ji tabelujeme s krokem Pi/2.

> k:=5:
    body:=[seq([a+i*(b-a)/(k-1),f(a+i*(b-a)/(k-1))],i=0..k-1)];

Z vybraných uzlových bodů vyrobíme příslušnou soustavu 5 rovnic o 5 neznámých ...

> A:=array(1..k,1..k):
    for i from 1 by 1 to k do
        A[i,1]:=1:
    end do:
    for i from 1 by 1 to k do
        for j from 2 by 1 to k do
            A[i,j]:=0:
        end do:
    end do:
    f(y):=1:
    for j from 2 by 1 to k do
        f(y):=f(y)*(y-body[j-1][1]):
        for i from j by 1 to k do
            A[i,j]:=eval(f(y),y=body[i][1]):
        end do:
    end do:
eval(A);

> B:=array(1..k):
    for i from 1 by 1 to k do
        B[i]:=body[i][2]:
    end do:
    eval(B);

... a vyřešíme.

> a:=linsolve(A,B);

Řešení opět dosadíme do vzorce ...

> N(x):=sum(a[n]*product((x-body[m][1]),m=1..n-1), n=1..k);

... a nakonec ještě zjdnodušíme.

> N(x):=expand(N(x));