CORRECTION: many thanks to Joel David Hamkins and Willie Wong; yes, we should take $I_1 = I_2 = 0$ (now edited out and removed). For some reason "constant" and "zero" keep getting mixed up together in my head!!
You can get lots of examples satisfying the vertical and horizontal integral conditions by choosing $f(x,y) = g(x)h(y)$ where $\int_0^1 g(x) dx = I_1$, $\int_0^1 \int_0^1 h(y) dy = I_2$ are arbitrary0$. The diagonal integrals being equal to $I_1 I_2$ 0$ gives you two equations relating $g,h$, but these still leave very many possibilities for $g,h$; then finallyg,h$.
(Note: I wrote $I_1$, $I_2$ before, but I've removed these to make it clearer).
Finally, taking arbitrary finite linear combinations of such $f$ (and also infinite linear combinations, if you are careful with convergence) gives you yet more examples.
If you want to restrict $f$ to, say, analytic functions, then this won't work - but you didn't say this in your question! (Although I suppose maybe you meant this, since you do say "analytic" magic square!)
EDIT: I'm definitely not claiming that every continuous example $f$ can be obtained in this way!! I just want to show that there are many, many possibilities for $f$.