xx : ([1.0, 2.0, 4.0, 6.0, 9.0]) $ 
yy : ([2.0, 4.0, -1.0, 3.0, 8.0]) $ 
n : length(xx) $
 xy : create_list([xx[i],yy[i]],i,1,n);
 l(i,x) := product((x-xx[j])/(xx[i]-xx[j]),j,1,i-1) *
           product((x-xx[j])/(xx[i]-xx[j]),j,i+1,n)$
 pol1(x) := sum(yy[i]*l(i,x),i,1,n);
pol1(x);
 expand(pol1(x));
in:pol1(2);
 wxplot2d( [[discrete,xy], pol1], [x, xx[1], xx[n]],
 [style, [points,3,2], [lines,1,1]]) $
 wxplot2d( [[discrete,xy],pol1], [x, xx[1], xx[n]], [y,-3,16],
[style, [points,3,1], [lines,1,2]]) $