%computer problem 4.2 #9
%in addition to printing the differences in values, 
%plots both the function and its interpolant

f = @(x) abs(x);
n = 11;

c = zeros(1,n);
y = zeros(1,n);

%first we interpolate f at 11 equally spaced nodes in [-1,1]

x = linspace(-1,1,n);      %equally spaced nodes in [-1,1]

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

for i=1:n
    c(i) = f(x(i));        %these values will be overwritten...
end

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

%now we compare our interpolant with exact values of f
%at 41 equally spaced points in [-1,1]

x_val = linspace(-1,1,41);
y_int = zeros(1,41);     %y_int = values of polynomial interpolant on x_val

for i=1:41
    y_int(i) = c(1);     %inner loop computes y_int(i) = P(x_val(i))
    p = 1;
    for j=2:n
        p = p*(x_val(i) - x(j-1));
        y_int(i) = y_int(i) + c(j)*p;
    end
end

y_val = zeros(1,41); 

for i=1:41
    y_val(i) = f(x_val(i));   %values of f on x_val
end

disp('difference of interpolating and exact values')
y_val - y_int

subplot(1,2,1)
plot(x_val,y_val,'-b',x_val,y_int,'-r',x,y,'ob')  %here we plot the difference
legend('exact','interpolant')
title('equally spaced nodes')


%now we repeat all of the above, but with Chebyshev nodes

for i=1:n
    x(i) = cos((2*i-1)*pi/(2*n));   %Chebyshev nodes on [-1,1]
end

for i=1:n
    y(i) = f(x(i));                 %y-values on Chebyshev nodes
end

for i=1:n
    c(i) = f(x(i));                 %these values will be overwritten...
end

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

%now we compare our interpolant with exact values of f
%at 41 equally spaced points in [-1,1]

x_val = linspace(-1,1,41);
y_int = zeros(1,41);     %y_int = values of polynomial interpolant on x_val

for i=1:41
    y_int(i) = c(1);     %inner loop computes y_int(i) = P(x_val(i))
    p = 1;
    for j=2:n
        p = p*(x_val(i) - x(j-1));
        y_int(i) = y_int(i) + c(j)*p;
    end
end

y_val = zeros(1,41); 

for i=1:41
    y_val(i) = f(x_val(i));   %values of f on x_val
end

disp('difference of interpolating and exact values')
y_val - y_int

subplot(1,2,2)   %plot the difference between exact and interpolant
plot(x_val,y_val,'-b',x_val,y_int,'-r',x,y,'ob')
legend('exact','interpolant')
title('chebyshev nodes')