
% Do a simple convergence study
N = 10;
h = []; E = [];
for j = 1:5
  h(j) = 1/(N+1);
  E(:,j) = particlebox(N);
  N = N*2;
end
loglog(h, pi^2-E(1,:));

% Rate of convergence
est_order = log( (pi^2-E(1,end-1))/(pi^2-E(1,end)) )/log(2);
fprintf('Estimated order of convergence: %f\n', est_order);
