% p10.m: Fourier collocation method for computing the whole stability
% spectrum of vortex solitons in 2D equation: iU_t+U_xx+U_yy+|U|^2U=0

 N=20; Lx=20; Ly=20; Nx=4*N; Ny=4*N; maxuerror=1e-7;
 dx=Lx/Nx;x=-Lx/2:dx:Lx/2-dx; kx0=2*pi/Lx;kx=[0:Nx/2-1 -Nx/2:-1]*kx0;
 dy=Ly/Ny;y=-Ly/2:dy:Ly/2-dy; ky0=2*pi/Ly;ky=[0:Ny/2-1 -Ny/2:-1]*ky0;
 [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=2; DT=0.95; uerror=1;
 while uerror > maxuerror          % 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;
 for j=1:2*N+1
   for k=1:2*N+1
     efactor=exp(-i*((j-N-1)*kx0*X+(k-N-1)*ky0*Y))/(Lx*Ly);
     c0(j,k)=dx*dy*sum(sum(G0.*efactor));
     c1(j,k)=dx*dy*sum(sum(G1.*efactor));
     c2(j,k)=dx*dy*sum(sum(G2.*efactor));
   end
 end
 D=zeros((2*N+1)^2); C0=D; C1=D; C2=D;
 for jj=1:(2*N+1)^2
   j=-N+floor((jj-1)/(2*N+1));  k=-N+mod(jj-1, 2*N+1);
   for kk=1:(2*N+1)^2
     m=-N+floor((kk-1)/(2*N+1));  n=-N+mod(kk-1, 2*N+1);
     ind1=j-m+N+1; ind2=k-n+N+1;
     if ind1 > 0 & ind1 < 2*N+2 & ind2 > 0 & ind2 < 2*N+2
        C0(jj, kk)=c0(ind1, ind2);
        C1(jj, kk)=c1(ind1, ind2);
        C2(jj, kk)=c2(ind1, ind2);
     end
   end
   D(jj, jj)=-((j*kx0)^2+(k*ky0)^2);
 end
 M=[ C0    D+C1
     D+C2   -C0 ];
 M=sparse(M); eigvalues=eigs(i*M, 200, 'SM');
 subplot(221)
 mesh(x, y, abs(u)); view([-30, 50]);
 subplot(222)
 plot(eigvalues, '.', 'markersize', 10); axis([-3 3 -3 3])