# Adaptive plotting of two-variable functions $z=f(x,y)$ algorithm pseudocode?

I am looking for explanations of algorithms to adaptively sample a function of two variables $f(x,y)$, in a given domain $x_0\le x \le x_1$, $y_0\le y \le y_1$. Intuitively, I want to sample more densely regions were $f(x,y)$ has more variation. I know almost nothing a priori about the function, except that inside the domain $x_0\le x \le x_1$, $y_0\le y \le y_1$, there maybe regions where it is infinite (that I obviously don't want to sample, so the algorithm should detect that too), and that everywhere else, the function is continuous. These singular regions are "smooth sets".