% compare sums with approximations
r=input('number>');

e1=[];
e2=[];

% Each data point is a relative error that becomes part of output vectors
% e1 and e2.  These vectors are plotted as a function of N.
for N=1:1:r
   x = 1:1:N;

   m = sum(x);
   p = N^2/2;
   e1 = [e1 (abs(m-p)/m)];

   m = sum(x.^2);
   p = N^3/3;
   e2 = [e2 (abs(m-p)/m)];
end
plot(1:r,e1,'b+',1:r,e2,'go');
title('Percentage error: sum((1:N)^i) vs. N^{(i+1)}/(i+1)');
ylabel('Error: abs(m-p)/p where m=sum(1:N)^i and p=N^{(i+1)}/(i+1)');
xlabel('N');
legend('sum(1:N)','sum((1:N)^2)');