
Lx=6*pi; Ly=Lx; N=64; max_iteration=5000; lambda_tol=1e-10;                 
dx=Lx/N; dy=Ly/N; x=-Lx/2:dx:(Lx/2-dx); y=-Ly/2:dy:(Ly/2-dy); 
kx=[0:N/2-1  -N/2:-1]*2*pi/Lx; ky=[0:N/2-1  -N/2:-1]*2*pi/Ly; 
[X,Y]=meshgrid(x,y); [KX,KY]=meshgrid(kx,ky); 
V=1.5*(cos(X).^2+cos(Y).^2); 
mu=0; U=1.7*(exp(-(X+pi/2).^2-(Y+pi/2).^2)+i*exp(-(X-pi/2).^2-(Y+pi/2).^2)-exp(-(X-pi/2).^2-(Y-pi/2).^2)-i*exp(-(X+pi/2).^2-(Y-pi/2).^2)); 
c=2; DTU=1; DT=0.4; 

for nn=1:max_iteration                % SOM for the vortex soliton
     Uold=U; 
     L0U=ifft2(-(KX.^2+KY.^2).*fft2(U))+(mu-V).*U+abs(U).^2.*U; 
     MinvL0U=ifft2(fft2(L0U)./(KX.^2+KY.^2+c)); 
     L1HermitMinvL0U=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0U))+(abs(U).^2+mu-V).*MinvL0U+2*U.*real(conj(U).*MinvL0U); 
     MinvL1HermitMinvL0U=ifft2(fft2(L1HermitMinvL0U)./(KX.^2+KY.^2+c));
     U=U-MinvL1HermitMinvL0U*DTU; 
     Uerror(nn)=sqrt(sum(sum(abs(U-Uold).^2))*dx*dy); 
     if Uerror(nn) < 1e-10
         break
     end
end

v=(exp(-(X+pi/2).^2-(Y+pi/2).^2)+exp(-(X-pi/2).^2-(Y-pi/2).^2))-(1+i)*(exp(-(X-pi/2).^2-(Y+pi/2).^2)+exp(-(X+pi/2).^2-(Y-pi/2).^2));
w=(exp(-(X-pi/2).^2-(Y+pi/2).^2)+exp(-(X+pi/2).^2-(Y-pi/2).^2))-(1+i)*(exp(-(X+pi/2).^2-(Y+pi/2).^2)+exp(-(X-pi/2).^2-(Y-pi/2).^2));

T12=mu+2*abs(U).^2-(U.^2+conj(U).^2)/2-V; 
T21=mu+2*abs(U).^2+(U.^2+conj(U).^2)/2-V; 
T11=(conj(U).^2-U.^2)/2; 
T22=(U.^2-conj(U).^2)/2; 

for nn=1:max_iteration                % MSOM for the eigenvalue
    vold=v; wold=w; 
    L10=-i*(T11.*v-ifft2((KX.^2+KY.^2).*fft2(w))+T12.*w); 
    L20=-i*(T22.*w-ifft2((KX.^2+KY.^2).*fft2(v))+T21.*v);
    Minvv=ifft2(fft2(v)./(c+KX.^2+KY.^2));  Minvw=ifft2(fft2(w)./(c+KX.^2+KY.^2));
    lambda=sum(sum(conj(Minvv).*L10+conj(Minvw).*L20))/sum(sum(conj(Minvv).*v+conj(Minvw).*w)); 
    Lambda(nn)=lambda; 
        
    L1=L10-lambda*v; 
    L2=L20-lambda*w;           
    MinvL1=ifft2(fft2(L1)./(c+KX.^2+KY.^2));   
    MinvL2=ifft2(fft2(L2)./(c+KX.^2+KY.^2));  

    LMinvL1=i*(conj(T11).*MinvL1-ifft2((KX.^2+KY.^2).*fft2(MinvL2))+conj(T21).*MinvL2)-conj(lambda)*MinvL1; 
    LMinvL2=i*(conj(T22).*MinvL2-ifft2((KX.^2+KY.^2).*fft2(MinvL1))+conj(T12).*MinvL1)-conj(lambda)*MinvL2; 
    MinvLMinvL1=ifft2(fft2(LMinvL1)./(c+KX.^2+KY.^2));   
    MinvLMinvL2=ifft2(fft2(LMinvL2)./(c+KX.^2+KY.^2));  
                    
    if nn == 1
        v=v-MinvLMinvL1*DT; 
        w=w-MinvLMinvL2*DT; 
    else
        L1G1=-i*(T11.*G1 -ifft2((KX.^2+KY.^2).*fft2(G2))+T12.*G2)-lambda*G1; 
        L1G2=-i*(T22.*G2 -ifft2((KX.^2+KY.^2).*fft2(G1))+T21.*G1)-lambda*G2;   
        MinvL1G1=ifft2(fft2(L1G1)./(c+KX.^2+KY.^2));
        MinvL1G2=ifft2(fft2(L1G2)./(c+KX.^2+KY.^2));
        MG1=ifft2(fft2(G1).*(c+KX.^2+KY.^2));
        MG2=ifft2(fft2(G2).*(c+KX.^2+KY.^2));
        alpha1=1/sum(sum(conj(MG1).*G1+conj(MG2).*G2)); 
        alpha2=1/(DT*sum(sum(conj(L1G1).*MinvL1G1+conj(L1G2).*MinvL1G2))); 
        alpha=alpha1-alpha2; 
        beta=sum(sum(conj(G1).*LMinvL1+conj(G2).*LMinvL2)); 
        v=v-(MinvLMinvL1-alpha*beta*G1)*DT; 
        w=w-(MinvLMinvL2-alpha*beta*G2)*DT; 
        lambda_error=abs(Lambda(nn)-Lambda(nn-1)); 
        format short e; [lambda  lambda_error]
        if lambda_error < lambda_tol
            break
        end
    end
    G1=v-vold; G2=w-wold; 
end
semilogy(abs(Lambda(2:end)-Lambda(1:end-1)),'linewidth', 2); 
xlabel('number of iterations'); ylabel('error')
