function ris=newton_NLS(m,F,x,power,index_max)
global mu sigma

beta=x(1);
cFix=x(2*index_max);
x(2*m+2:end)=[];
Xk=x(2:end);

Mat=[];
for k=1:m
    Matk=[];
    for l=1:m
        Matk=[Matk,[-2,1;1,-2]*F(k,l)];
    end
    Mat=[Mat;Matk];
end
Mat=sigma*Mat+diag(reshape(repmat((pi*[1:1:m]).^2,[2,1]),[],1));
MatL=Mat+diag(repmat([sigma*mu-beta,sigma*mu+beta],[1,m]));
fm=MatL*Xk;
%Newton scheme to determine the numerical solution of the eigenvalue/eigenvector system
x(2*index_max)=[];

attempts=0;
while norm(fm,inf)>10^(-power)
    Dfm=[((-1).^[1:1:2*m])'.*Xk,MatL(:,1:2*(index_max-1)),MatL(:,2*index_max:end)];   % non derivo rispetto c2
  
    x1=x-Dfm\fm;
    x=x1;
    beta=x(1);
  
    Xk=[x(2:2*index_max-1);cFix;x(2*index_max:end)];
    MatL=Mat+diag(repmat([sigma*mu-beta,sigma*mu+beta],[1,m]));
    fm=MatL*Xk;    
   attempts=attempts+1;
   if attempts==20
       disp('The newton procedure for the eigenvalues does not converge')
       return
   end
end
ris=[x(1:2*index_max-1);cFix;x(2*index_max:end)];

end