ibox=initboxes2d(n,z)

function ibox=initboxes2d(n,z)
% Initial coarse discretization on z x z
% where z = [a b], a < b

% Initial interior centers
deltax = (z(2)-z(1))/n;
x = linspace(z(1),z(2),n+1)';
x = x(1:n) + 0.5*deltax;
[X,Y] = meshgrid(x,x);
xy = num2cell([X(:) Y(:)],2);
[ibox(1:n^2).midpt] = deal(xy{:});
[ibox.dimen] = deal([deltax deltax]);
[ibox.isboundary] = deal(false);
[ibox.child] = deal([]);
[ibox.Nchild] = deal(0);
[ibox.depth] = deal(0);

% South centers
xy = num2cell([x z(1)*ones(n,1)],2);
[sbox(1:n).midpt] = deal(xy{:});
[sbox.dimen] = deal([deltax 0]);
[sbox.isboundary] = deal(true);
[sbox.child] = deal([]);
[sbox.Nchild] = deal(0);
[sbox.depth] = deal(0);

% North centers
xy = num2cell([x z(2)*ones(n,1)],2);
[nbox(1:n).midpt] = deal(xy{:});
[nbox.dimen] = deal([deltax 0]);
[nbox.isboundary] = deal(true);
[nbox.child] = deal([]);
[nbox.Nchild] = deal(0);
[nbox.depth] = deal(0);

% East centers
xy = num2cell([z(2)*ones(n,1) x],2);
[ebox(1:n).midpt] = deal(xy{:});
[ebox.dimen] = deal([0 deltax]);
[ebox.isboundary] = deal(true);
[ebox.child] = deal([]);
[ebox.Nchild] = deal(0);
[ebox.depth] = deal(0);

% West centers
xy = num2cell([z(1)*ones(n,1) x],2);
[wbox(1:n).midpt] = deal(xy{:});
[wbox.dimen] = deal([0 deltax]);
[wbox.isboundary] = deal(true);
[wbox.child] = deal([]);
[wbox.Nchild] = deal(0);
[wbox.depth] = deal(0);

ibox = [ibox sbox nbox ebox wbox];