f = @(x,y) (3*x^2-y^3-x);
g = @(x,y) (2*y^2-2*x-y);
x0=1.2; y0=3.4;

Jac = @(x,y) [6*x-1, -3*y^2;
                 -2, 4*y-1;];

J=Jac(x0,y0);
disp("J(x0,y0)=");
disp(J);

disp("det(J(x0,y0))=");
disp(det(J));

J_Inv=inv(J);
disp("J_inverse(x0,y0)=");
disp(J_Inv);

disp("( f(x0,y0) , g(x0,y0) )=");
disp([f(x0,y0); g(x0,y0)]);

ans=[x0; y0;]-(J_Inv*[f(x0,y0); g(x0,y0)]);
disp("(x1,y1)=");
disp(ans);
