def METHOD_NAME(p, dx, dy, b):
    pn = np.empty_like(p)
    pn = p.copy()
    for q in range(nit):
        pn = p.copy()
        p[1:-1, 1:-1] = ((pn[2:, 1:-1]+pn[0:-2, 1:-1])*dy**2 +
                         (pn[1:-1, 2:]+pn[1:-1, 0:-2])*dx**2) / \
                        (2*(dx**2+dy**2)) - \
            dx**2*dy**2/(2*(dx**2+dy**2))*b[1:-1, 1:-1]
        p[-1, :] = p[-2, :]   # dp/dy = 0 at y = 2
        p[0, :] = p[1, :]   # dp/dy = 0 at y = 0
        p[:, 0] = p[:, 1]   # dp/dx = 0 at x = 0
        p[:, -1] = 0        # p = 0 at x = 2
    return p