x=-5:(0.01):5;
y=1./(1+3*x.^2);
hold on
plot(x,y,'b-');
xx=(-5):(0.5):5;  
yy=1./(1+3*xx.^2); 
x=(-5):(0.01):5; 
n=length(x);   
y=zeros(n,1);  
for ii=1:21
    y1=ones(n,1);
    for jj=1:21
        if  ii ~= jj 
            y1=y1.*(x-xx(jj))'/(xx(ii)-xx(jj));
        end
    end
    y=y+y1*yy(ii);
end
hold on              
plot(x,y,'g-');  
xx=(-5):(0.5):5;
yy=1./(1+3*xx.^2);
hold on
plot(xx,yy,'r*');
x=-5:(0.5):5;
y=1./(1+3*x.^2);
xx=-5:(0.01):5;
yy=spline(x,y,xx);
hold on
plot(xx,yy,'k');