function r = ConvOrder(Error)
% This function makes a plot of log(error(n+1)) vs log(error(n)),
% given an array of error values, and calculates the order of
% convergence using the first and last data points.

n = prod(size(Error));           % determine size of array "Error"
if  (Error(n) == 0)  n = n - 1;  % throw out last point if error==0
end;

logError = log(Error(1:n));
X = logError(1:(n-1));
Y = logError(2:n);
order = round( (Y(n-1) - Y(1)) / (X(n-1) - X(1)) );
C = exp( (sum(Y) - order * sum(X)) / (n-1) );  % Optional
Line = log(C) + order * X;
tt = sprintf('Convergence order = %d; C=%.3g\n', order, C);

plot(X, Y, 'ob');
hold on;
set(gca, 'fontsize', 16); 
plot(X, Line, '-r');
title(tt);
xlabel('ln( \epsilon_n )'); ylabel('ln( \epsilon_{n+1} )'); 
legend('ln(\epsilon_{n+1}) vs. ln(\epsilon_n)', 'Line fit', 'location', 'best');
