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.$