function [x,it]=newtonit(F,x0,t,tol)
xold=x0;
H=jacobian(F,xold,t);
invH=inv(H);
xnew=xold-inv(H)*F(xold,t);
it=1;
while norm(xnew-xold)>tol
  it=it+1;
  xold=xnew;  
  H=jacobian(F,xold,t);
  xnew=xold-inv(H)*F(xold,t);
end
x=xnew;
endfunction
