It depends on your initial guesses. Some work, some don't (not unusual for Newton's method!): Also, better practice to set D to be the Jacobian, rather than its inverse, then use backslash division in the iteration see below:
x = 2; y = 2;
[x,y] = newton3(x,y);
disp([x y])
disp([x^4+y^4-2*x*y^5 x^6+x^2+y^4-4 ])
function [x, y] = newton3(x,y)
for N = 1:30
D = [4*x^3-2*y^5 4*y^3-10*x*y^4; 6*x^5+2*x 4*y^3]; %%%%%%%
f = x^4+y^4-2*x*y^5 ;
g = x^6+x^2+y^4-4 ;
z = [x y]' ;
z = z - D\[f g]' ; %%%%%%%%
x = z(1);
y = z(2);
end
end
