TPlane := proc(f,x,y,x0,y0)
local fx,fy,z0;
fx := subs(x=x0,y=y0, diff(f,x)):
fy := subs(x=x0,y=y0, diff(f,y)):
z0 := subs(x=x0,y=y0, f):
# the tangent plane is then given by z =

sort(fx*(x-x0) + fy*(y-y0) + z0,[x,y]);
end: