x = linspace(1,30,10);   // x values
y = rand(x);  // y values
 
function[P]=lagrange(X,Y) 
    P=0;
    for i=1:length(X)
        L=1;
        for j=[1:i-1,i+1:length(X)] 
            L=L*(poly(0,"x")-X(j))/(X(i)-X(j));
        end
        P=P+L*Y(i); 
    end
endfunction

clf()       // clear plot window
n = 100     // number of points in interpolated curve
xx = linspace(1,30,n); // get n evenly spaced values in interval
yspline = interp(xx, x, y, splin(x,y,"natural"));  // interpolate at all x values (builds vector)
 
P = lagrange(x,y);
ylagrange = horner(P,xx);
 
plot2d(xx, ylagrange, 3 )
plot2d(xx, yspline,2 )  // plot curve throught interpolated points
plot2d(x,y,-9)   
legends(["Cubic Spline","Lagrange","Random Data Point"], [2 3 -9], "ur", %f)
xtitle("Cubic Spline vs Lagrange Interpolation")
show_window()
