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 Andrey's solution for continuous functions will still work 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.$