I wondered whether it holds, for a function $f:A\subset\mathbb{R}^3\to\mathbb{R}$, $f\in C(A)$, that, for all $(x_0,y_0,z_0)\in A$, $$\lim_{\sqrt{h_x^2+h_y^2+h_z^2}\to 0}\frac{1}{h_x h_y h_z}\int_{z_0}^{z_0+h_z}\int_{y_0}^{y_0+h_y}\int_{x_0}^{x_0+h_x}f(x,y,z)dxdydz=f(x_0,y_0,z_0)$$I think that, if such an identity held, it would allow us to define quantities like the densities used in physics.
Obviously, for all $\varepsilon>0$ there exists a $\delta$ such that, if $\sqrt{h_x^2+h_y^2+h_z^2}<\delta$ (and therefore $|h_x|,|h_y|,|h_z|$ $<\delta$), then $$\Bigg|\frac{1}{h_z}\int_{z_0}^{z_0+h_z}\int_{y_0}^{y_0+h_y}\int_{x_0}^{x_0+h_x}f(x,y,z)dxdydz-\int_{y_0}^{y_0+h_y}\int_{x_0}^{x_0+h_x}f(x,y,z_0)dxdy\Bigg|<\varepsilon$$$$\Bigg|\frac{1}{h_y}\int_{y_0}^{y_0+h_y}\int_{x_0}^{x_0+h_x}f(x,y,z_0)dxdy-\int_{x_0}^{x_0+h_x}f(x,y_0,z_0)dx\Bigg|<\varepsilon$$$$\Bigg|\frac{1}{h_x}\int_{x_0}^{x_0+h_x}f(x,y,z_0)dxdydz-f(x_0,y_0,z_0)\Bigg|<\varepsilon$$but I am not able to use these inequalities to find an arbitrarily small $\tilde{\varepsilon}$ to majorate $\big|\frac{1}{h_xh_yh_z}\int_{z_0}^{z_0+h_z}\int_{y_0}^{y_0+h_y}\int_{x_0}^{x_0+h_x}f(x,y,z)dxdydz-f(x_0,y_0,z_0)\big|$ with.
Does what I am trying to prove hold and, if it does, how can it be proved? I $\infty$-ly thank anybody for any answer.
You don't need Lebesgue density / differentiation theorem, as $f$ is continuous.
As I am lazy, I will prove it in 2D and assume without loss of generality $x_0 = y_0 = 0$ and $f(0,0) = 0$. We need to show that $$ \lim_{r\to 0} \sup_{\|h\| < r, h_xh_y\ne 0} \left | \frac1{h_x h_y} \int_0^{h_y} \int_0^{h_x} f(x,y) \; \mathrm dx\mathrm dy \right | = 0.$$
Now, let $\epsilon > 0$. Then, by continuity of $f$ there is some $r > 0$ such that for every $h$ with $\|h\|< r$ it follows $|f(h_x,h_y)| < \epsilon.$ Notice that for every $(x,y)$ with $|x| \le |h_x|$ and $|y|\le |h_y|$ we have $$ \sqrt{x^2 + y^2} \le \| h \| < r. $$ Thus, for $h_xh_y \ne 0$ and $s_x = \operatorname{sign}(h_x), s_y = \operatorname{sign}(h_y)$ we obtain \begin{align} \left | \frac1{h_x h_y} \int_0^{h_y} \int_0^{h_x} f(x,y) \; \mathrm dx\mathrm dy \right | &= \frac1{|h_x| |h_y|} \left | \int_0^{|h_y|} \int_0^{|h_x|} f( s_x x, s_y y) \; \mathrm dx\mathrm dy \right | \\ &\le \frac1{|h_x| |h_y|} \int_0^{|h_y|} \int_0^{|h_x|} \left | f( s_x x, s_y y)\right | \; \mathrm dx\mathrm dy \\ &< \frac1{|h_x| |h_y|} \int_0^{|h_y|} \int_0^{|h_x|} \epsilon \; \mathrm dx\mathrm dy \\ &= \epsilon. \end{align}