function seq = BR_Newton_MP2(tol,a0,b0,max_iter)
% Newton method in N-dim
% seq = BR_Newton_MP2(0.00000001,0.5,0.5,10000);

a = a0;
b = b0;
f = BR_f_MP(a,b);
seq = [a0,b0,f,0];
H = zeros(2,2);
DEL = zeros(2,1);
iter = 0;
while tol<f&(iter<max_iter)
    iter = iter + 1;
    H(1,1) = 4 - 4*b + 2*b^2;
    H(1,2) = 0.6 - 4*a +4*a*b;
    H(2,1) = H(2,1);
    H(2,2) = 2*a^2;
    DEL(1) = 4*a + 0.6*b - 4*a*b + 2*a*b^2 - 2;
    DEL(2) = 0.6*a - 2*a^2 + 2*(a^2)*b;
    del_theta = atan2(DEL(2),DEL(1));
    s = H\-DEL;
    a = a + s(1);
    b = b + s(2);
    fs = BR_f_MP(a,b);
    f = fs;
    seq = [seq; a,b,f,180*del_theta/pi];
end
