function [X,Y,iter] = NEWTONmultiD(eps,T,A,F,Niter)
% [X,Y,iter] = NEWTONmultiD(eps,T,A,F,Niter); 
JF = jacobian(F,T);  
X = A; e  = 2*eps; iter = 0; 
while iter < Niter && e > eps  
    U = eval(subs(F,T,X)); V = eval(subs(JF,T,X)); 
    if abs(det(V)) < eps
        X = X; e = 0;  
    else 
        Xnew = X' - V^(-1)*U; e = norm(Xnew-X'); X = Xnew'; 
    end
iter = iter +1;     
end
Y = U';
end 


