Plot on a 2D grid a spherically symmetric function
No need to iterate if the function works with arrays, and you use broadcasting to populate a grid:

def fn(X,Y):
    r = np.sqrt(X**2+Y**2)
    return 1/r**2
fn(grid_x[:,None], grid_y[None,:])


np.meshgrid and np.mgrid (or np.ogrid) can also be used to form the X,Y grid.  

fn(*np.meshgrid(grid_x, grid_y))
fn(*np.ogrid[-Grid_side:Grid_side:N_pix*1j, -Grid_side:Grid_side:N_pix*1j])


Look up their docs, and practice.  Practice with simple things like grid_x[:,None]+grid_y as well.