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

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

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

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

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

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