f1=x^3 + y^3 - 2 x y;
f2=x^2 + y^3-1;

ContourPlot[{f1==0,f2==0},{x,-3,3},{y,-3,3}]

Solve[{f1==0,f2==0},{x,y},Reals]//N

jac=D[{f1,f2},{{x,y}}]

ContourPlot[{f1==0,f2==0,Det[jac]==0},{x,-3,3},{y,-3,3}]

newton[{x0_,y0_}]:=({x,y}-Inverse[jac].{f1,f2})/.{x->x0,y->y0}

newton[{a,b}]//Simplify

NestList[newton,{1.,1.},9]//TableForm

newton2[{x_,y_}]:=N[newton[{x,y}],30];

NestList[newton2,{1,1},13]//TableForm