iVar = input('Enter an r to compute sum(k^r) for k=1 to n: r =');
r=iVar;
% Since we have proven that the series:
%  sum(k^i,k=1..n) is an (i+1)^st degree polynomial, we can use linear
%  algebra to solve for the polynomial coefficients. Let:
%  f(x) = sum(k^i,k=1..x) for x an integer. Since f(x) is an (i+1)^st
%  degree polynomial, it can be uniquely defined by i+1 distinct points.
%  Thus we calculate f(1), f(2), ..., f(i+1). We have:
%   f(1) = a(1)*1^(i+1) + a(2)*1^i + ... + a(i)*1^1 + a(i+1)
%   f(2) = a(1)*2^(i+1) + a(2)*2^i + ... + a(i)*2^1 + a(i+1)
%     .  .  .                          
%   f(i+1) = a(1)*(i+1)^(i+1) + a(2)*(i+1)^i + ... + a(i)*(i+1)^1 + a(i+1)
%  This forms a linear system which we can then solve.
%
%  Example: For i=1 we have a quadratic. We compute:
%   f(1) = a(1)+a(2)+a(3) = 1
%   f(2) = a(1)*2^2+a(2)2^1+a(3) = 1+2 = 3
%  f(3) = a(1)*3^2+a(2)*3^1+a(3) = 1+2+3 = 6
 
% The matrix we need to compute has a name, it is the Vandermonde matrix.
% Thus we can just use the function for a Vandermonde matrix, named vander
j=1:iVar+2;
 
A = vander(j);
A
b= cumsum(j.^iVar)'; % This means "cumulative sum"
 
 
coeffs = A\b;
 
% Display: This part is just to display all of the polynomials "nicely"
fprintf('The polynomial is: ');
j=1;
if(abs(coeffs(j)) >= 1e-3)
    fprintf('%5.4gx^%d',coeffs(j),iVar+2-j);
end
for j=2:iVar
    if(abs(coeffs(j)) >= 1e-3 && coeffs(j) > 0 )
        fprintf('+%5.4gx^%d',coeffs(j),iVar+2-j);
    elseif(abs(coeffs(j)) >= 1e-3 && coeffs(j) < 0 )
        fprintf('%5.4gx^%d',coeffs(j),iVar+2-j);
    end
end
j=iVar+1;
if(abs(coeffs(j)) >= 1e-3 && coeffs(j) > 0 )
    fprintf('+%5.4gx',coeffs(j));
elseif(abs(coeffs(j)) >= 1e-3 && coeffs(j) < 0 )
    fprintf('%5.4gx',coeffs(j));
end
 
j=iVar+2;
if(abs(coeffs(j)) >= 1e-3 && coeffs(j) > 0 )
    fprintf('+%5.4g',coeffs(j));
elseif(abs(coeffs(j)) >= 1e-3 && coeffs(j) < 0 )
    fprintf('%5.4g',coeffs(j));
end
fprintf('\n');