This post is my attempt to answer the question from here .

The procedure **ContoursWithLabels** has 2 required parameters: **Expr** is an expression in **x** and **y** variables, **Range1** and **Range2** are ranges for **x** and **y** . In this case, the output is the list of floats for the contours and 8 black contours (with labels) (the axis of coordinates as a box).

The optional parameters: **Number** is positive integer - the number of contours (by default Number=8), **S** is a set of real numbers C for contours (for which Expr=C) (by default S={}), **GraphicOptions** is a list of graphic options for plotting (by default GraphicOptions=[color = black, axes = box]), **Coloring** is an equality Coloring=list of color options for **plots[dencityplot]** command (by default Coloring=NULL).

The code of the procedure:

**restart;**

**ContoursWithLabels := proc (Expr, Range1::(range(realcons)), Range2::(range(realcons)), Number::posint := 8, S::(set(realcons)) := {}, GraphicOptions::list := [color = black, axes = box], Coloring::`=` := NULL)**

**local r1, r2, L, f, L1, h, S1, P, P1, r, M, C, T, p, p1, m, n, A, B, E;**

**uses plots, plottools;**

**f := unapply(Expr, x, y);**

**if S = {} then r1 := rand(convert(Range1, float)); r2 := rand(convert(Range2, float));**

**L := [seq([r1(), r2()], i = 1 .. 205)];**

**L1 := convert(sort(select(a->type(a, realcons), [seq(f(op(t)), t = L)]), (a, b) ->is(abs(a) < abs(b))), set);**

**h := (L1[-6]-L1[1])/Number;**

**S1 := [seq(L1[1]+(1/2)*h+h*(n-1), n = 1 .. Number)] else**

**S1 := convert(S, list) fi;**

**print(Contours = evalf[2](S1));**

**r := k->rand(20 .. k-20); M := []; T := [];**

**for C in S1 do**

**P := implicitplot(Expr = C, x = Range1, y = Range2, op(GraphicOptions), gridrefine = 3);**

**P1 := [getdata(P)];**

**for p in P1 do**

**p1 := convert(p[3], listlist); n := nops(p1);**

**if n < 500 then m := `if`(40 < n, (r(n))(), round((1/2)*n)); M := `if`(40 < n, [op(M), p1[1 .. m-11], p1[m+11 .. n]], [op(M), p1]); T := [op(T), [op(p1[m]), evalf[2](C)]] else**

**if 500 <= n then h := floor((1/2)*n); m := (r(h))(); M := [op(M), p1[1 .. m-11], p1[m+11 .. m+h-11], p1[m+h+11 .. n]]; T := [op(T), [op(p1[m]), evalf[2](C)], [op(p1[m+h]), evalf[2](C)]]**

**fi; fi; od; od;**

**A := plot(M, op(GraphicOptions));**

**B := plots:-textplot(T);**

**if Coloring = NULL then E := NULL else E := ([plots:-densityplot])(Expr, x = Range1, y = Range2, op(rhs(Coloring))) fi;**

**display(E, A, B);**

**end proc:**

** **

Examples of use:

**ContoursWithLabels(x^2+y^2, -3 .. 3, -3 .. 3);**

** **

**ContoursWithLabels(x^2-y^2, -5 .. 5, -5 .. 5, {-20, -15, -10, -5, 0, 5, 10, 15, 20}, [color = black, thickness = 2, axes = box], Coloring = [colorstyle = HUE, colorscheme = ["White", "Red"], style = surface]);**

** **

The next example, I took from here:

**ContoursWithLabels(sin(1.3*x)*cos(.9*y)+cos(.8*x)*sin(1.9*y)+cos(.2*x*y), -5 .. 0, 2 .. 5, {seq(-2 .. 2, 0.5)}, [color = black, axes = box], Coloring = [colorstyle = HUE, colorscheme = ["Cyan", "Red"], style = surface]);**

** **

There are many more examples can be found in the attached file.

ContoursWithLabels1.mw

Edit. The attached file has been corrected.