Suppose $f: \mathbb{R}\times\mathbb{R}\to\mathbb{R}$ is has continuous partial derivatives and
$$4f(x,y)=f(x+\delta,y+\delta)+f(x-\delta,y+\delta)+f(x-\delta,y-\delta) + f(x+\delta,y-\delta)$$
for all $(x,y)$ in $\mathbb{R}\times\mathbb{R}$ and all $\delta$ in $\mathbb{R}$.
I don't believe that $f$ is necessarly harmonic but I cannot construct a counter-example. Is $f$ harmonic?
Yes, the Taylor series works. Actually $C^2$ suffices for the remainder term, although my sophomore calculus book gives the proof using $C^3.$ I get $$ 4 f(x_0, y_0) = 4 f(x_0, y_0) + \left( 2 f_{xx}(x_0, y_0) + 2 f_{yy}(x_0, y_0) \right) \delta^2 \; + \; o( \delta^2 ) $$ and $$ \left( 2 f_{xx}(x_0, y_0) + 2 f_{yy}(x_0, y_0) \right) \delta^2 \; = \; o( \delta^2 ) $$ and $$ 2 \left( f_{xx}(x_0, y_0) + f_{yy}(x_0, y_0) \right) \; = \; 0 $$
LATER EDIT: unless I am vastly mistaken this argument still works if we put in the caveat $ | \delta | < \Delta = \Delta(x_0, y_0), $ that is we only require your equation for small $\delta$ and even say that the allowable size of $\delta$ depends on the position of the center point that I am calling $(x_0, y_0).$ But with this change we can build an easy discontinuous example of your relation, take $$ f(x_0, y_0) = 1, \; \; if \; \; y_0 > 0, $$ $$ f(x_0, y_0) = 0, \; \; if \; \; y_0 = 0, $$ $$ f(x_0, y_0) = -1, \; \; if \; \; y_0 < 0. $$
Then your relation holds for $ | \delta | < | y_0 | $ when $y_0 \neq 0$ and holds for all $\delta$ when $ y_0 = 0.$
$f$ is harmonic under the weaker assumption that it is just continuous.
Multiplying the identity $$f(x+\delta,y+\delta)+f(x-\delta,y+\delta)+u(x-\delta,y-\delta) + f(x+\delta,y-\delta)-4f(x,y)=0$$ with a test function $g\in C_0^{2}(\mathbb R^2)$ and integrating the result over $\mathbb R^2$, it's easy to see that $$\int_{\mathbb R^2}\left(g(x+\delta,y+\delta)+g(x-\delta,y+\delta)+g(x-\delta,y-\delta) + g(x+\delta,y-\delta)-4g(x,y)\right)f(x,y)\ dxdy=0.$$ Applying to the latter equality the argument in Will's answer, we obtain that, for every $g\in C_0^{2}(\mathbb R^2)$ $$\int_{\mathbb R^2}(g_{xx}+g_{yy})f\ dxdy=0. \qquad\qquad (*)$$
By a theorem of Kellog, every continuous solution $f$ to $(*)$ satisfies the mean value property $$f(P)=\frac{1}{2\pi}\int_{0}^{2\pi}f(P+re^{i\phi})d\phi$$ for any $P\in \mathbb R^2$ and any $r>0$. Therefore, $f$ is harmonic.
Update. As Pietro Majer and BS indicated, the argument works for locally integrable $f$. This and the comment of Will Jagy above make me think that any discontinuous function that solves the discrete Laplace equation for every $\delta>0$ should have highly pathological properties. Perhaps, non-measurable solutions to Cauchy's functional equation $$f(x+y)=f(x)+f(y)$$ might be a close analogy.
To reduce the case of a continuous $f:\mathbb{R}^2\to\mathbb{R}$ to, say the $C^2$ or $C^3$ case, one can simply mollify $f$ via convolution with a smooth kernel with compact support. Then $f_\epsilon:=f*\phi_\epsilon$ still enjoys the "N-S-W-E mean property" above, so it's harmonic as already seen in previous answers, and since as $\epsilon\to0$ the $f_\epsilon$ converge to $f$ uniformly on compact sets together with all second order derivatives, $f$ is harmonic too. Note that the analogous holds for a continuous $f$ (or even just locally integrable, it works as well) on an open subset of $\mathbb{R}^n$.
Variation: take $\phi$ as above and moreover with radial symmetry. Then $f*\phi$ is harmonic as before, so $f*\phi*\phi=f*\phi$ because harmonic functions are invariant by convolution with radial symmetric kernels (it's just a weighted radially symmetric mean value property). Since $\phi$ has compact support we can simplify $\phi$ in the last equality (this is standard via Fourier transform) and get $f*\phi=f$ so f itself is harmonic.
