function [x,nit,flag]=newtonsistemi(F,JF,x,p)
nmax=100;
toll=1.e-12;
flag=1;
step=0;
JFx=feval(JF,x);
[L,U,P]=lu(JFx);
for nit=1:nmax
    step=step+1;
    if step==p
        step=0;
        JFx=feval(JF,x);
        [L,U,P]=lu(JFx);
    end
    Fx=feval(F,x);
    z1=L\(P*Fx);
    z=U\z1;
    x=x-z;
    if norm(z)<toll, flag=0; return, end
end
