Let $u(x,y)$ be a smooth function on the square $S:=[0, 1] \times [0, 1]$ (see, for example, Wiki for the definitions) and $\varepsilon > 0$.
Is it possible to approximate $u(x,y)$ by a function $h(x,y)$ harmonic on $S$ s.t. $$\max_{(x,y)\in S} |u(x,y)-h(x,y)|\le \varepsilon ? $$
No, this is not possible because (locally) uniform limits of harmonic functions are harmonic so the limit $u=\lim\limits_{\epsilon\to 0} h_\epsilon$ would have to be harmonic to start with.
Well, I guess it depends on what one means by "harmonic", but at least in the flat Euclidean setting there is no amibguity. In this basic setting the fact that uniform limits of harmonic functions are harmonic immediately follows from the standard characterization of harmonic functions by the mean-value property (which is trivially stable under uniform limits).
A counterexample is therefore given by any non-harmonic smooth function. (I will not insult MO's readership by giving an explicit such function!)
Think of it like this: In dimension 1 harmonic functions are affine, and clearly it is impossible to approximate an arbitrary smooth function uniformly (or in any reasonable topology, for that matters) by sequences of affine functions. Going to higher dimensions does not help with the "approximability".
