xf=6*pi; yf=xf; N=64; max_iteration=5000; lambda_tol=1e-10; 
x=-xf/2:xf/N:(xf/2-xf/N); y=-yf/2:yf/N:(yf/2-yf/N); dx=xf/N; dy=yf/N; 
kx=[0:N/2-1  -N/2:-1]*2*pi/xf; ky=[0:N/2-1  -N/2:-1]*2*pi/yf; 
[X,Y]=meshgrid(x,y); [KX,KY]=meshgrid(kx,ky); dkx=2*pi/xf; dky=2*pi/yf;
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)); 
DTU=1; DTlambda=2; c=2; lambda=0; 

for nn=1:max_iteration                % SOM for the vortex soliton
     Uold=U; 
     L0U=ifft2(-(KX.^2+KY.^2).*fft2(U))-V.*U+abs(U).^2.*U+mu*U; 
     MinvL0U=ifft2(fft2(L0U)./(KX.^2+KY.^2+c)); 
     LMinvL0U=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0U))+(abs(U).^2+mu-V).*MinvL0U+2*U.*real(conj(U).*MinvL0U); 
     MinvLMinvL0U=ifft2(fft2(LMinvL0U)./(KX.^2+KY.^2+c));
     U=U-MinvLMinvL0U*DTU; 
     Uerror=sqrt(trapz(trapz(abs(U-Uold).^2))*dx*dy); 
     if Uerror < 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                % OOM for the eigenvalue
        vold=v; wold=w; lambdaold=lambda; 
        Lv=-i*(T11.*v - ifft2((KX.^2+KY.^2).*fft2(w))+T12.*w); 
        Lw=-i*(T22.*w - ifft2((KX.^2+KY.^2).*fft2(v))+T21.*v);
        lambda=trapz(trapz(Lv.*conj(v)+Lw.*conj(w)))/trapz(trapz(abs(v).^2+abs(w).^2)); 
        Lambda(nn)=lambda; 
        L1=Lv-lambda*v; 
        L2=Lw-lambda*w; 
        dv1=ifft2(fft2(Lv-lambda*v)./(c+KX.^2+KY.^2));  Tv=v+0.5*dv1*DTlambda; 
        dw1=ifft2(fft2(Lw-lambda*w)./(c+KX.^2+KY.^2));  Tw=w+0.5*dw1*DTlambda; 

        Lv=-i*(T11.*Tv - ifft2((KX.^2+KY.^2).*fft2(Tw))+T12.*Tw); 
        Lw=-i*(T22.*Tw - ifft2((KX.^2+KY.^2).*fft2(Tv))+T21.*Tv);   
        lambda=trapz(trapz(Lv.*conj(Tv)+Lw.*conj(Tw)))/trapz(trapz(abs(Tv).^2+abs(Tw).^2)); 
        dv2=ifft2(fft2(Lv-lambda*Tv)./(c+KX.^2+KY.^2));  Tv=v+0.5*dv2*DTlambda; 
        dw2=ifft2(fft2(Lw-lambda*Tw)./(c+KX.^2+KY.^2));  Tw=w+0.5*dw2*DTlambda; 

        Lv=-i*(T11.*Tv - ifft2((KX.^2+KY.^2).*fft2(Tw))+T12.*Tw); 
        Lw=-i*(T22.*Tw - ifft2((KX.^2+KY.^2).*fft2(Tv))+T21.*Tv);   
        lambda=trapz(trapz(Lv.*conj(Tv)+Lw.*conj(Tw)))/trapz(trapz(abs(Tv).^2+abs(Tw).^2)); 
        dv3=ifft2(fft2(Lv-lambda*Tv)./(c+KX.^2+KY.^2));  Tv=v+dv3*DTlambda; 
        dw3=ifft2(fft2(Lw-lambda*Tw)./(c+KX.^2+KY.^2));  Tw=w+dw3*DTlambda; 

        Lv=-i*(T11.*Tv - ifft2((KX.^2+KY.^2).*fft2(Tw))+T12.*Tw); 
        Lw=-i*(T22.*Tw - ifft2((KX.^2+KY.^2).*fft2(Tv))+T21.*Tv);   
        lambda=trapz(trapz(Lv.*conj(Tv)+Lw.*conj(Tw)))/trapz(trapz(abs(Tv).^2+abs(Tw).^2)); 
        dv4=ifft2(fft2(Lv-lambda*Tv)./(c+KX.^2+KY.^2));  v=v+(dv1+2*dv2+2*dv3+dv4)*DTlambda/6; 
        dw4=ifft2(fft2(Lw-lambda*Tw)./(c+KX.^2+KY.^2));  w=w+(dw1+2*dw2+2*dw3+dw4)*DTlambda/6; 
        
        lambda_error=abs(lambda-lambdaold); 
        format short e; [lambda  lambda_error]
        if lambda_error < lambda_tol
            break
        end
end
semilogy(abs(Lambda(2:end)-Lambda(1:end-1)), 'r--', 'linewidth', 2); 
xlabel('number of iterations'); ylabel('error')