C..1 2次元のNewton法

f[x_,y_]:={x^2-y^2+x+1,2 x y +y}
Df[a_,b_]:=Module[
           {x,y},
           Transpose[{D[f[x,y],x],D[f[x,y],y]}] /. {x->a,y->b}
          ]
{x,y}={1,1}
Do[{x,y}={x,y}-Inverse[Df[x,y]].f[x,y];
   Print[{x,y},"=",N[{p,q},20]],
   {6}
  ]

あるいは、変数の方もベクトル的にリストを使って書いて、

f[{x_,y_}]:={x^2-y^2+x+1,2 x y +y}
Df[{a_,b_}]:=Module[
           {x,y},
           Transpose[{D[f[{x,y}],x],D[f[{x,y}],y]}] /. {x->a,y->b}
          ]
xk={1,1}
Do[xk=xk-Inverse[Df[xk]].f[xk]; Print[xk,"=",N[xk,20]], {6}]