%interpolates Runge function f at 5,9,13, and 17 Chebyshev node
%points in [-5,5]; interpolation polynomial is plotted against f

f = @(x) (1+x^2)^(-1);

for m=1:4
    
n = 4*m;
x = zeros(1,n+1);
for i=1:n+1
    x(i) = 5*cos((2*i-1)*pi/(2*n+2));   %Chebyshev nodes on [-5,5]
end
c = zeros(1,n+1);              %coefficients of Newton interpolant

for i=1:n+1
    c(i) = f(x(i));
end

for j=2:n+1
    for i=n+1:-1:j
        c(i) = (c(i)-c(i-1))/(x(i)-x(i-j+1));
    end
end

x_new = linspace(-5,5,1000);
z = zeros(1,1000);   %z = values of polynomial interpolant on x_new

for i=1:1000
    z(i) = c(1);     %inner loop computes z(i) = P(x_new(i))
    p = 1;
    for j=2:n+1
        p = p*(x_new(i) - x(j-1));
        z(i) = z(i) + c(j)*p;
    end
end

y_new = zeros(1,1000); 

for i=1:1000
    y_new(i) = f(x_new(i));   %values of f on x_new
end

y = zeros(1,n+1);
for i=1:n+1
    y(i) = f(x(i));
end

if m==1
    subplot(2,2,1)
end
if m==2
    subplot(2,2,2)
end
if m==3
    subplot(2,2,3)
end
if m==4
    subplot(2,2,4)
end

plot(x_new,y_new,'-b',x_new,z,'-r',x,y,'o')

end

