syms x x1 y n m s Q
fplot('1/(1+x^2)',[-5,5])
hold on;
n=3;
x=linspace(-5,5,50);
for m=1:3;
    Q=0;
    s=1;
    n=n*5;
    x1=linspace(-5,5,n);
    y=1./(1+x1.^2);
    for k=1:50
    for j=1:n;
        for i=1:n;
            if i~=j;
              s(k)=((x(k)-x1(i))*s(k))/(x1(j)-x1(i))*y(i);
          end
      end
      Q(k)=Q(k)+s(k);
  end
end
  plot(x,Q)
  hold on;
end

hold on;