Here's a simple vector valued function:
In: f[{x_, y_}] := {x^2 y + 1, E^x + y}
The Jacobian:
In: jacobian[f, {x, y}]
Out: {{2 x y, x^2}, {E^x, 1}}
To construct a Jacobian function that we can evaluate we have to create a Mathematica. To evaluate the derivatives in place we use other vector component names and replace them after calculating the derivative:
In: df[{x_, y_}] := jacobian[f, {a, b}] /. {a -> x, b -> y}
With the Jacobian function it's easy to set up the iteration with FixedPointList and an anonymous function:
In: FixedPointList[(# - Inverse[df[#]].f[#]) &, {-0.7, -2.5}, 20]
Out: {{-0.7, -2.5}, ... , {0.703467, -2.02075}}
-a
