Create jacobian Function from a vector Function

Question

This is related to this question and maybe this question. I wanted a function that takes an n-parameter vector function and outputs a Jacobian matrix function, also n-parameter. I think this works:

J[f_]:=Module[
   {x=Unique[],v},
   v=Thread[Subscript[x, Range[Length[{##}]]]];
   Outer[D,f@@v,v]/.Thread[v->{##}]
]&

J[Function[{x,y},{x y+x^2-2,Cos[x-y]Sin[x]-y}]][1.5,2.5]

D[{x y+x^2-2,Cos[x-y]Sin[x]-y},{{x,y}}]/.Thread[{x,y}->{1.5,2.5}]

but I was wondering if there is a cleaner way that doesn't use rules. I'm posting this question even though it's basically solved, because it took me a while, as do most problems involving manipulating pure functions.

Answer 1

J[f_] := Transpose[Map[Function[j, j[##]],Map[Apply[Derivative[##][f] &, #] &, IdentityMatrix[Length[{##}]]]]] &

J[{#1 #2 + #1^2 - 2, Cos[#1 - #2] Sin[#1] - #2} &][1.5, 2.5]

(* {{5.5, 1.5}, {0.877583, -1.83936}} *)

Or in more readable form

J[f_] := Module[{n, j},
  n = Length[{##}];
  j = Map[Apply[Derivative[##][f] &, #] &, IdentityMatrix[n]];
  Transpose[Map[Function[j, j[##]], j]]
] &

... but I will only claim to have met one of your requirements. Incidentally, do you also observe that there's no way to have a named slot sequence in a pure function?