r := (x,y) -> x*i+y*j+f(x,y)*k; 'fx(a,b)' = subs(x=a,diff(r(x,b),x)); 'fy(a,b)' = subs(y=b,diff(r(a,y),y)); N := crossprod([1,0,D1f(a,b)],[0,1,D2f(a,b)]); TPlane := innerprod(N,evalm([x,y,z]-[a,b,f(a,b)])) = 0; readlib(isolate): isolate(TPlane,z); collect(",[D1f(a,b),D2f(a,b)]); sort(",f(a,b));