format long e

p0=[-2 3]';
f1=@(x)x(1)+2*x(2)+0.03*sin(x(1)+x(2))-4;
f2=@(x)5*x(1)+6*x(2)+0.07*cos(x(1)-x(2))-8;
f=@(x)[f1(x) f2(x)]';
df1x1=@(x)1+0.03*cos(x(1)+x(2));
df1x2=@(x)2+0.03*cos(x(1)+x(2));
df2x1=@(x)5-0.07*sin(x(1)-x(2));
df2x2=@(x)6+0.07*sin(x(1)-x(2));
J=@(x)[df1x1(x) df1x2(x);df2x1(x) df2x2(x)];

for i=1:1000
  p=p0-J(p0)\f(p0);
  if abs(p-p0)<10^-15
    fprintf("x1=%.15f x2=%.15f with %d interation\n",p0(1),p0(2),i);
    break
    end
    p0=p;
    end