function f = newtonmv_f(v)
x = v(1);
y = v(2);
f = zeros(2, 1);
f(1) = x^2 + y^2 - 4;
f(2) = x*y - 1;
% plot solution
xi = -2:.02:2;
xi1 = -2.2:.02:-0.1;
xi2 = 0.1:.02:2.2;
plot(xi,sqrt(4-xi.^2),'r',xi,-sqrt(4-xi.^2),'r',...
    xi1,1./xi1,'b',xi2,1./xi2,'b')
axis([-2.2 2.2 -2.2 2.2])
