Some example plots

plot3d(sin(x+y)*(x-y), [x,-5,5], [y,-5,5], [plot_format, gnuplot])$

here the points outside the domain are given wrong

plot3d(sqrt(2-x^2-y^2), [x,-2,2], [y,-2,2], [plot_format, gnuplot])$

there we go:

draw3d(xlabel="x", ylabel="y", zlabel="z",
       title="sqrt(2-x^2-y^2)",
       line_width=2,
       parametric_surface(v*cos(t),v*sin(t),
                          sqrt(2-v^2),
                          v,0,sqrt(2),t,0,2*%pi),
       surface_hide=true,
       terminal=wxt)$

Here's a sliding trace (for the slider, select
the graph by clicking on it and then move the slider in the
toolbar)

with_slider_draw3d(
       a, makelist(i/2, i, -10, 10),
       surface_hide=false,
       xlabel="x", ylabel="y", zlabel="z",
       line_width=1,
       explicit(sin(x+y)*(x-y), x,-5,5, y,-5,5),
       color=blue,
       nticks=50,
       line_width=4,
       pic_width=700,
       pic_height=500,
       parametric(t,a,sin(t+a)*(t-a),t,-5,5),
       line_width=2,
       color=green,
       points_joined=true,
       point_size=0,
       points([[-5,a,-9],[5,a,-9],[5,a,9],[-5,a,9],[-5,a,-9]]),
       color=black,
       label([concat("y=",string(a)),5,a,-9]),
       user_preamble="set xyplane at -11")$

Plot with contours, and then just the contours.  First a function with
a maximum, then a saddle

draw3d(xlabel="x", ylabel="y", zlabel="z",
       title="z=exp(-(sin(x)^2+y^2))",
       line_width=2,
       explicit(exp(-(sin(x)^2+y^2)), x,-2.5,2.5, y,-2.5,2.5),
       contour_levels=10,
       contour=both,
       surface_hide=true,
       terminal=wxt)$

draw3d(xlabel="x", ylabel="y", zlabel="z",
       title="z=exp(-(sin(x)^2+y^2))",       
       line_width=2,
       explicit(exp(-(sin(x)^2+y^2)), x,-2.5,2.5, y,-2.5,2.5),
       contour_levels=10,
       contour=map,
       surface_hide=true,
       terminal=wxt)$

draw3d(xlabel="x", ylabel="y", zlabel="z",
       title="z=x^2-y^2",
       line_width=2,
       explicit(x^2-y^2, x,-2.5,2.5, y,-2.5,2.5),
       contour_levels=10,
       contour=both,
       surface_hide=true,
       terminal=wxt)$

draw3d(xlabel="x", ylabel="y", zlabel="z",
       title="z=x^2-y^2",
       line_width=2,
       explicit(x^2-y^2, x,-2.5,2.5, y,-2.5,2.5),
       contour_levels=10,
       contour=map,
       surface_hide=true,
       terminal=wxt)$

Level sets of a 3 dimensional function (for the slider, select
the graph by clicking on it and then move the slider in the
toolbar)

with_slider_draw3d(
       a, makelist(i, i, -15, 10),
       surface_hide=true,
       xlabel="x", ylabel="y", zlabel="z",
       enhanced3d=true,
       line_width=2,
       pic_width=700,
       pic_height=500,
       implicit(-(z*2)^2+x^2-y^2=a, x,-5,5, y,-5,5, z,-5,5),
       title=concat("x^2-y^2-(2z)^2=",string(a)),
       user_preamble= "set xyplane at -5")$