clear 
a = -1;
b = 1;
n = 12;

f = @(x) 1./(1+25.*x.^2);
%f = @(x) sin(2*pi*x);

x1 = linspace(a,b,n+1);
y1 = f(x1);
p1 = polyfit(x1,y1,n);

x2 = -cos(((0:n) + 0.5)*pi/(n+1));
y2 = f(x2);
p2 = polyfit(x2,y2,n);

x_fine = linspace(a,b,300);
yp1 = polyval(p1,x_fine);
yp2 = polyval(p2,x_fine);

plot(x1,y1,'o',x2,y2,'x',x_fine,yp1,x_fine,yp2);