% the Cartesian grid
GRID.x = linspace(ax,bx,nx);
GRID.y = linspace(ay,by,ny);

GRID.hx = GRID.x(2)-GRID.x(1); GRID.hy = GRID.y(2)-GRID.y(1);
GRID.n = nx;
GRID.nx = nx;
GRID.ny = ny;

[GRID.X,GRID.Y] = ndgrid(GRID.x,GRID.y);

% linear index of boundary points (of the computational domain)
GRID.bdry = [[ones(ny,1)    [1:ny]'];...
	     [nx*ones(ny,1) [1:ny]'];...
	     [[1:nx]' ones(nx,1)];...
	     [[1:nx]' ny*ones(nx,1)]];
GRID.bdry = sub2ind([nx,ny],GRID.bdry(:,1),GRID.bdry(:,2));
GRID.bdry = unique(GRID.bdry);

clear nx ny hx hy