OK is this insanity?

%Code written by Echo Liu

format 'long';

k = 10;
NMax = 50000;
N = 10;
dummy = 0;
dummy1 = 0;
dummy2 = 0;
%QhatN = zeros((NMax-k)+1,1);
%a = zeros(N,1);
IiVec = zeros(NMax,1);
Qn = zeros(NMax,1);
a = zeros(NMax+1,1);
QhatN = zeros(NMax-2,1);
InteEval = 0;
QhatDiff = zeros(NMax-3,1);
QnDiff = zeros(NMax-1,1);
%f = @(r) r^(-1)*(cos(r));
f = @(r) r.^-1.*(sin(r.^-1.*log(r)));
a(1) = 1;
xaxisforQhatDiff = 1:1:NMax-3;
xaxisforQnDiff = 1:1:NMax-1;
xaxisforQn = 1:1:NMax;
xaxisforQhatN = 1:1:NMax-2;
for i=1:NMax
  b = fzero(@(x) x*exp(x)-i*pi,0);
  a(i+1) = exp(-b);

  dummy = a(i);
  dummy1 = a(i+1);
  dummy2 =a(i) - a(i+1);

 InteEval = integral(f,a(i),a(i+1));
 IiVec(i) = InteEval;
   %if (i == 1)
    %   Qn(i) = IiVec(i);
   %else
    %   Qn(i) = Qn(i-1) + IiVec(i);
   %end

end

%Summing up all the elements
Qn(1) = IiVec(1);
for s=2:NMax
    Qn(s) = Qn(s-1) + IiVec(s);
end

%Aitken Part
for r=1:NMax-2
    %QhatN(r) = QhatN(r)- (((IiVec(r+1))^(2))/(IiVec(r+2))-(IiVec(r+1)));
    QhatN(r) = Qn(r) - (((Qn(r+1)-Qn(r))^(2))/(Qn(r+2)-2*Qn(r+1)+Qn(r)));
end

for z=1:length(QhatN)-1
    QhatDiff(z) = QhatN(z+1) - QhatN(z);
end

for f=1:length(Qn)-1
    QnDiff(f) = Qn(f+1)-Qn(f);
end
% we take the abs value because we
% are interested in the abs error for each p(i)

figure;
plot(log10(xaxisforQnDiff),log10(abs(QnDiff')))
title(['Graph of differences up to N for regular version for N = ' num2str(NMax) '.'])
xlabel(['log of N'])
ylabel(['log of differences'])

% we take the abs value because we
% are interested in the abs error for each p(i)
figure;
plot(log10(xaxisforQhatDiff),log(abs(QhatDiff')))
title(['Graph of differences up to N for Aitkens method for N = ' num2str(NMax) ','])
xlabel(['log of N'])
ylabel(['log of differences'])

figure;
plot(xaxisforQn,Qn');
title(['Qn vs N']);

figure;
plot(xaxisforQhatN,QhatN');
title(['QhatN vs N']);

figure;
plot(log10(xaxisforQnDiff),log10(abs(QnDiff')),'r',log10(xaxisforQhatDiff),log(abs(QhatDiff')),'b')
%title{'Differences for both sequences over N with log base 10 both axes')