%naive algorithm for interpolating f(x) = (1-x^N)/(1-x)
%fails for large N due to numerical errors

N = 10;
f = @(x) (1-x^N)/(1-x);

%first we set up the linear system to solve for coefficients
%of the interpolating polynomial for f at nodes x(i), i=1,...,N

x = zeros(1,N);
for i=1:N
    x(i) = 1+i;
end

A = zeros(N,N+1);
for i=1:N
    for j=1:N
        A(i,j) = x(i)^(j-1);   %A = Vandermonde matrix for nodes x(i)
    end
end

for i=1:N
    A(i,N+1) = f(x(i));        %augmented Vandermonde matrix
end

[m,n] = size(A);       %A is m by n

%now we solve the linear system corresponding to A via naive Gaussian elimination

%forward elimination

for k=1:m-1
    for i=k+1:m
        c = A(i,k)/A(k,k);   %c = multiple of row i to subtract from row j
        for j=1:n
            A(i,j) = A(i,j)-c*A(k,j);
        end
    end
end

%backward elimination

for k=m:-1:2          %loop goes from m to 2 (the -1 means it's backwards)
    for i=1:k-1
        c = A(i,k)/A(k,k);   %c = multiple of row i to subtract from row j
        for j=1:n
            A(i,j) = A(i,j)-c*A(k,j);
        end
    end
end

%put into reduced row echelon form

for i=1:m
    c = A(i,i);       %c = matrix entry at pivot in row i
    for j=1:n
        A(i,j) = A(i,j)/c;    %divide row i by c
    end
end

disp('the reduced matrix is...')
A

%the far right column of A should be the all ones vector; why?
