clear 
format long

df1dx1 = @(x) 1 + 0.03*cos(x(1)+x(2));
df1dx2 = @(x) 2 + 0.03*cos(x(1)+x(2));
df2dx1 = @(x) 5 - 0.07*sin(x(1)-x(2));
df2dx2 = @(x) 6 + 0.07*sin(x(1)-x(2));

J = @(x) [df1dx1(x) df1dx2(x) ; df2dx1(x) df2dx2(x)]; 

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;

G = @(x) x - J(x)\[f1(x) f2(x)]';

p0 = [-2,3]';


for i=1:50
  p = p0;
  p0 = G(p);
end 

f1(p0)
f2(p0)

p0