function [s]=find_node(x,y,xmax,xmin,nx,ymin,ymax,ny,G);
    hx=(xmax-xmin)/(nx-1);
    hy=(ymax-ymin)/(ny-1);
    s=find(G.xy(:,1)>=x - hx & G.xy(:,1)<=x & ...
           G.xy(:,2)>=y - hy & G.xy(:,2)<=y);
end