% p11.m:  Newton-conjugate-gradient method for computing eigenvalues
% of vortex solitons in the 2D NLS equation: iU_t+U_xx+U_yy+|U|^2U=0

 Lx=30; Ly=30; Nx=128; Ny=128; maxerror=1e-10; errorCG=1e-2;
 x=-Lx/2:Lx/Nx:Lx/2-Lx/Nx; kx=[0:Nx/2-1  -Nx/2:-1]*2*pi/Lx;
 y=-Ly/2:Ly/Ny:Ly/2-Ly/Ny; ky=[0:Ny/2-1  -Ny/2:-1]*2*pi/Ly;
 [X, Y]=meshgrid(x, y); [TH R]=cart2pol(X,Y);
 [KX, KY]=meshgrid(kx, ky); K2=KX.^2+KY.^2;
 mu=1; u=2.5*R.*sech(R).*exp(i*TH); c=4; DT=1.95; uerror=1;
 while uerror > maxerror          % SOM for vortex soliton starts
    F=u.*conj(u)-mu;
    L0u=ifft2(-K2.*fft2(u))+F.*u;
    uerror=max(max(abs(L0u)))
    MinvL0u=ifft2(fft2(L0u)./(c+K2));
    L1MinvL0u=ifft2(-K2.*fft2(MinvL0u))+F.*MinvL0u   ...
              +2*u.*real(conj(u).*MinvL0u);
    MinvL1MinvL0u=ifft2(fft2(L1MinvL0u)./(c+K2));
    u=u-MinvL1MinvL0u*DT;
 end                               % SOM for vortex soliton ends
 G0=(u.*u-conj(u).*conj(u))/2;
 G1=-mu+2*u.*conj(u)-(u.*u+conj(u).*conj(u))/2;
 G2=-mu+2*u.*conj(u)+(u.*u+conj(u).*conj(u))/2;
 f1=tanh(0.5*R).*sech((0.5-0.2i)*R).^2;
 f2=tanh(0.5*R).*sech((0.5+0.5i)*R).^2;
 v=f1.*exp(-i*TH)+f2.*exp(-3*i*TH);           % initial conditions
 w=f1.*exp(-i*TH)-f2.*exp(-3*i*TH);
 ncg=0; nnt=0; fftM=(c+K2).^2;
 while 1                  % Newton-CG method for eigenvalue starts
    nnt=nnt+1;
    v2w2=sum(sum(conj(v).*v+conj(w).*w));
    Lv=i*( G0.*v-ifft2(K2.*fft2(w))+G1.*w);
    Lw=i*(-G0.*w-ifft2(K2.*fft2(v))+G2.*v);
    lambda=sum(sum(conj(v).*Lv+conj(w).*Lw))/v2w2;
    curlL0v=Lv-lambda*v;
    curlL0w=Lw-lambda*w;
    vwerror(nnt)=max(max(abs(curlL0v)+abs(curlL0w))); vwerror(nnt)
    numcg(nnt)=ncg;
    if vwerror(nnt) < maxerror
        break
    end
    L1HcurlL0v=-i*( conj(G0).*curlL0v-ifft2(K2.*fft2(curlL0w)) ...
                   +conj(G2).*curlL0w)-conj(lambda)*curlL0v;
    L1HcurlL0w=-i*(-conj(G0).*curlL0w-ifft2(K2.*fft2(curlL0v)) ...
                   +conj(G1).*curlL0v)-conj(lambda)*curlL0w;
    dv=0*v;  dw=0*w;                       % CG iterations start
    Rv=-L1HcurlL0v;  Rw=-L1HcurlL0w;
    MinvRv=ifft2(fft2(Rv)./fftM);  MinvRw=ifft2(fft2(Rw)./fftM);
    R2=sum(sum(conj(Rv).*MinvRv+conj(Rw).*MinvRw));  R20=R2;
    Dv=MinvRv;  Dw=MinvRw;
    while (R2 > R20*errorCG^2)
        L1Dv=i*( G0.*Dv-ifft2(K2.*fft2(Dw))+G1.*Dw)-lambda*Dv;
        L1Dw=i*(-G0.*Dw-ifft2(K2.*fft2(Dv))+G2.*Dv)-lambda*Dw;
        gamma=sum(sum(conj(v).*L1Dv+conj(w).*L1Dw))/v2w2;
        curlL1Dv=L1Dv-gamma*v;
        curlL1Dw=L1Dw-gamma*w;
        L1HcurlL1Dv=-i*( conj(G0).*curlL1Dv-ifft2(K2.*     ...
          fft2(curlL1Dw))+conj(G2).*curlL1Dw)-conj(lambda)*curlL1Dv;
        L1HcurlL1Dw=-i*(-conj(G0).*curlL1Dw-ifft2(K2.*     ...
          fft2(curlL1Dv))+conj(G1).*curlL1Dv)-conj(lambda)*curlL1Dw;
        a=R2/sum(sum(conj(Dv).*L1HcurlL1Dv+conj(Dw).*L1HcurlL1Dw));
        dv=dv+a*Dv;  dw=dw+a*Dw;
        Rv=Rv-a*L1HcurlL1Dv;  Rw=Rw-a*L1HcurlL1Dw;
        MinvRv=ifft2(fft2(Rv)./fftM);  MinvRw=ifft2(fft2(Rw)./fftM);
        R2old=R2;
        R2=sum(sum(conj(Rv).*MinvRv+conj(Rw).*MinvRw));
        b=R2/R2old;
        Dv=MinvRv+b*Dv;  Dw=MinvRw+b*Dw;
        ncg=ncg+1;
    end                                    % CG iterations end
    v=v+dv;  w=w+dw;
 end                      % Newton-CG method for eigenvalue ends
 lambda
 subplot(221); imagesc(x,y,abs(u)); axis xy
 subplot(222); semilogy(numcg, vwerror, 'o', numcg, vwerror);
 xlabel('number of CG iterations'); ylabel('error');
 axis([0 90 1e-12 100]);
 subplot(223); imagesc(x, y, abs(v)); axis xy
 subplot(224); imagesc(x, y, angle(v)); axis xy