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.

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.