% Plots the function and the polynomial interpolation

close all; clear all;

z = -5:.01:5;
plot(z,1./(1+z.^2)); % plot the function
hold on;

% define the table for the (x,y) values
l=0:20;
x = 5 * cos((2*l+1)/(42)*pi);
% x = -5:0.5:5
y = 1./(1+x.^2)

plot(x,y,'ro') % plot the table using red circles

m = size(x,2); % find the size of the table
a=x(1);
b=x(m);
n=200; % sample the polynomial interpolation at n-many points

% main iteration

for i=1:n+1  % loop for sampling points
  
  xi(i) = a+(i-1)*(b-a)/n; 
  yi(i) = 0;
  
  for j = 1:m; % loop for evaluation poly. interp. at x=xi(i) 
        
    prod = 1;
    for k=1:m
      if ( k ~= j)
	prod = prod * ( xi(i) - x(k) ) / ( x(j) - x(k) );
      end
    end
    
    yi(i) = yi(i) + y(j) * prod;
    
  end
    
end

plot(xi,yi); % plotting the ploynomial interpolation 
hold off;

% x = 5*cos( (2*(0:10)+1) / (22) * pi );  % Chebyshev nodes