There are many published algorithms.
Are you working from a function you can evaluate at any point, or
a set of data points which you must interpolate? Working from a
continuous, differentiable function isn't too hard. Turning
randomly spaced data points into clean contours is not well
defined.
One simple approach is to subdivide. Impose a grid on the
data and look at the points in each cell. If the Z range within
a cell is below a contour Z unit, that cell is more or less flat.
If not, subdivide the cell into 4 cells and recurse. Now
you have a quadtree of flat cells, each with a Z value.
John Nagle