# Create jacobian Function from a vector Function

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},
]&

J[Function[{x,y},{x y+x^2-2,Cos[x-y]Sin[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.

-

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}} *)


J[f_] := Module[{n, j},