Newtons metod flervariabelanalys

Jag försöker använda Newtons metod för att lösa ett system av ekvationerJag har kommit så långt som att rita ut noll-nivåkurvorna och har därmed hittat två approximativa lösningar.

Jag har kört fast vid att försöka hitta bättre lösningar mha Newtons metod, min lösning hamnar långt ifrån den riktiga lösningen, och jag förstår inte varför... Här är min kod och plot:

%% Uppgift 1
clc
clear
% Börjar med att rita upp noll-nivåkurvorna
% till f1 och f2, där dem skär varandra finns 
% lösningewn till f(x)=0

% Beskriv komponentfunktionerna f1 och f2.

% Låter dem bero på två separata variabler
% för att kunna beräkna dem med matriser av x1-
% och x2-värden.

% f(x) = (x1^3+x2^2-1 , exp(x1*x2)+x1+x2-2) 
% =>

f1=@(x1,x2)x1.^3+x2^2-1;
f2=@(x1,x2)exp(x1.*x2)+x1+x2-2;

% countour med vektorn [0 0], 
% ger vi 0 får vi ingen nivå

x1=linspace(-20,10,30); x2=linspace(-20,10,30);
[X1,X2]=meshgrid(x1,x2);
Z1=f1(X1,X2); Z2=f2(X1,X2);
contour(x1,x2,Z1,[0 0], 'g')
hold on
contour(x1,x2,Z2,[0 0], 'b')
grid on

% ser att f1 skär f2 vid två punkter, 
% x=(-3.5,-0.5)
% x=(0.5,-0.5)

% 1 - Newton

%(x1^3+x2^2-1 , exp(x1*x2)+x1+x2-2) 
f=@(x)[x(1).^3+x(2).^2-1 
    exp(x(1).*x(2)+x(1)+x(2)-2)];

Df=@(x)[3.*x(1).^2 2.*x(2)
    exp(x(1)+x(2))+1 exp(x(1)+x(2))+1];

x=[-3.45;-0.34];

kmax=10; tol=0.5e-8;
for k=1:kmax
    h=-Df(x)\f(x);
    x=x+h;
    disp([x' norm(h)])
    if norm(h)<tol, break, end
end

plot(x(1),x(2),'ro')