z := solve(A*x+B*y+C*z+D,z); SqrDist := (x-x0)^2+(y-y0)^2+(z-z0)^2; solve({diff(SqrDist,x), diff(SqrDist,y)},{x,y}); f := (x,y) -> x^2 + 2*x*y + 3*y^2: gf := grad(f(x,y),[x,y]); solve({gf[1],gf[2]},{x,y});