Consider a vector field v$(x,y)=xy$i +$(x^2+y^2)$j
- Now I tried to do the following , but I cannot see why this limit does not come equal to $3y_1$ ?? Please help.. I would be very thankful for your precious time! Hope it benefits you too ! :)
It was just a calculation mistake.. but the question now :-
Can we do this for any function with domain $R^2$ that outputs a 2-D vector, i.e. , integrate then take limit to find the divergence? Is there a proof in general case?

There are some mistakes in your derivation. You oriented the square as it was a circulation problem and it is not. Else you have to define the outwards vector perpendicular to the boundary line to calculate the flux. For the sides parallel to $x$ the vector is $\mathbf j$ for the upper line and $-\mathbf j$ for the lower (so, the flux, e. g. through the lower side is $q(x,y_1-a)\mathbf j·(-\mathbf j)=-q(x,y_1-a)$). Similarly $\mathbf i$ and $-\mathbf i$. Furher, the area of the quadrilateral you defined is $4ab$ (the variables run from $-a$ to $a$ and from $-b$ to $b$, giving sides of length $2a$ and $2b$). The integrals for the flux per unit of area are:
$$I(a,b, x_1,y_1)=$$
$$=\left(\dfrac{1}{4ab}\right)\left(\int_{x_1-b}^{x_1+b}\left(q(x,y_1+a)-q(x,y_1-a)\right)dx+\int_{y_1-a}^{y_1+a}\left(p(x_1+b,y)-p(x_1-b,y)\right)dy\right)=$$
$$=\left(\dfrac{1}{4ab}\right)\left(\int_{x_1-b}^{x_1+b}\left(x^2+(y_1+a)^2-x^2-(y_1-a)^2\right)dx+\int_{y_1-a}^{y_1+a}\left((x_1+b)y-(x_1-b)y\right)dy\right)=$$
$$=\left(\dfrac{1}{4ab}\right)\left(\int_{x_1-b}^{x_1+b}4y_1a\,dx+\int_{y_1-a}^{y_1+a}2bydy\right)=$$
$$=\left(\dfrac{1}{4ab}\right)\left(\big[4y_1ax\big]_{x_1-b}^{x_1+b}+\left[by^2\right]_{y_1-a}^{y_1+a}\right)=$$
$$=\left(\dfrac{1}{4ab}\right)\left(8y_1ab+4y_1ab\right)=3y_1$$
Not depending on $a$ or $b$, the limit yields the same.