Laplacian as limit of Integral Identity Let $\psi(\vec{r})$ be a scalar field,
show that:
$$\nabla^2 \psi(\vec{r})=\lim_{\rho \to 0} \frac{3}{\pi \rho^2} \int_\Omega \psi(\vec{r}')-\psi(\vec{r})d\Omega'$$
where $\rho=|\vec{r}-\vec{r'}|$, $\Omega$ is a sphere of radius $\rho$ with it's center at $\vec{r}$ (and $d\Omega'$ symbolises the solid angle)
what I have tried it a Taylor expansion of $\psi(\vec{r}')$ around $\vec{r}$
that is:
$$\psi(\vec{r}')=\psi(\vec{r})+\nabla\psi(\vec{r})\cdot (\vec{r}-\vec{r}')+..$$
which gives me (because the limith should kill all the higher terms in the expansion, I think..)
$$\nabla^2 \psi(\vec{r})=\lim_{\rho \to 0} \frac{3}{\pi \rho^2} \int_\Omega \nabla\psi(\vec{r})\cdot\vec{\rho} d\Omega'$$
now I know that divergence is defined as
$$\nabla \cdot \vec{F}=\lim_{\rho \to 0} \frac{1}{4\pi \rho^2} \int_\Omega \vec{F}\cdot \vec{dA}$$
so I get:
$$\lim_{\rho \to 0} \frac{1}{4\pi \rho^2} \int_\Omega\nabla\psi(\vec{r})\cdot \vec{dA}=\lim_{\rho \to 0} \frac{3}{\pi \rho^2} \int_\Omega \nabla\psi(\vec{r})\cdot\vec{\rho} d\Omega'$$
I just don't know how to show that these two limits are the same.
can someone help (if what I have done so far is correct of course)
 A: Firstly, there are several similar posts regarding that kind of problem: for example, this, this and this.
Secondly, there are some mistakes in your reasoning I think. Here is what stuff is.
Let us begin from calculating the following expression:
$$
I := \iint\limits_{S^{2}(\vec{r}, \rho)}\big(\psi(\vec{r}') - \psi(\vec{r})\big)\mathrm{d}\Omega,\label{eq:star}\tag{*}
$$ 
where $S^{2}(\vec{r}, \rho)$ denotes 2D-sphere with radius $\rho$ centered at $\vec{r}$, $\Omega$ - solid angle (and not a surface element so it doesn't include $\rho$).
Let us expand $\psi(\vec{r}')$ around $\vec{r}$:
$$
\psi(\vec{r}') = \psi(\vec{r}) + (\vec{r}' - \vec{r})\cdot\nabla\psi\big|_{\vec{r}} + \frac12\big((\vec{r}' - \vec{r})\cdot\nabla\big|_{\vec{r}}\big)^2\psi + o(\rho^2)\label{eq:expand}\tag{1}
$$
NOTE: that in this formula, for example, vector $\nabla\psi$ is taken at the point $\vec{r}$ and is constant vector - it doesn't depend on $\vec{r}'$. Substitution of \eqref{eq:expand} into \eqref{eq:star} gives
$$
\iint\limits_{S^{2}(\vec{r}, \rho)}\big((\vec{r}' - \vec{r})\cdot\nabla\psi\big|_{\vec{r}} + \frac12\big((\vec{r}' - \vec{r})\cdot\nabla\big|_{\vec{r}}\big)^2\psi + o(\rho^2)\big)\mathrm{d}\Omega
$$
The first term disappears because of its antisymmetry along point $\vec{r}$. Let's investigate the second term. It could be shown (not a hard exercise) that
$$
\big((\vec{r}' - \vec{r})\cdot\nabla\big|_{\vec{r}}\big)^2\psi = \sum_{k=1}^{3}(r'_k - r_k)^2\frac{\partial^2\psi}{\partial x_k^2}\bigg|_{\vec{r}} + \sum_{i\neq j}^3(r'_i - r_i)(r'_j - r_j)\frac{\partial^2\psi}{\partial x_i\partial x_j}\bigg|_{\vec{r}}.\label{eq:psi2}\tag{2}
$$
Integrating this term (with $\frac12$ multiplier) gives $\frac{2\pi\rho^2}{3}$. NOTE: the second sum in \eqref{eq:psi2} due to its antisymmetry is also disappears when integrated.
The third term in \eqref{eq:star} is $o(\rho^2)$ when integrated so we have 
$$
\Delta\psi(\vec{r}) = \lim_{\rho\to 0} \frac{3}{2\pi\rho^2}I.
$$ 
