function [glx,gly,gx,gy] = grid_2d (ax,bx,ay,by,Nx,Ny)

%-------------------
% Generate a 2D grid
%-------------------

%-----
% prepare
%-----

Dx = (bx-ax)/Nx;
Dy = (by-ay)/Ny;

%---
% generate the grid lines
%---

 for i=1:Nx+1
   glx(i) = ax+(i-1.0)*Dx;
   plot([glx(i),glx(i)],[ay,by],'c')
 end

 for j=1:Ny+1
   gly(j) = ay+(j-1.0)*Dy;
   plot([ax,bx],[gly(j),gly(j)],'c')
 end

%---
% generate the grid points
%---

 for i=1:Nx+1
   for j=1:Ny+1
     gx(i,j) = glx(i);
     gy(i,j) = gly(j);
   end
 end

%---
% done
%---

return