
%Note: For Eqs. (142)-(143), 
%  L1 = [dxx+dyy+V-mu             U
%           0.5U          0.5dxx+2.5dyy-0.5-mu]
%  The Hermitian of L1 is the transpose of L1


Lx=50; Ly=50; N=64;                         % grid parameters
max_iteration=1e4; error_tolerance=1e-10; 
x=-Lx/2:Lx/N:Lx/2-Lx/N; dx=Lx/N; kx=[0:N/2-1  -N/2:-1]*2*pi/Lx;
y=-Ly/2:Ly/N:Ly/2-Ly/N; dy=Ly/N; ky=[0:N/2-1  -N/2:-1]*2*pi/Ly;
[X, Y]=meshgrid(x, y); [KX, KY]=meshgrid(kx, ky);  

mu=0.1; DT=0.59; 
U=sech(0.35*sqrt(X.^2+Y.^2));               % initial conditions
V=0.3*sech(0.5*sqrt(X.^2+Y.^2)); 
      
for nn=1:max_iteration                      % MSOM iterations start
    Uold=U; 
    Vold=V; 
    L0U=ifft2(-(KX.^2+KY.^2).*fft2(U))+U.*V-mu*U; 
    L0V=0.5*(ifft2(-(KX.^2+5*KY.^2).*fft2(V))+0.5*U.^2-V)-mu*V; 
    MinvL0U=ifft2(fft2(L0U)./(mu+KX.^2+KY.^2));  
    MinvL0V=ifft2(fft2(L0V)./(mu+0.5+0.5*KX.^2+2.5*KY.^2));           
    L1HermitMinvL0U=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0U))+(V-mu).*MinvL0U+0.5*U.*MinvL0V; 
    L1HermitMinvL0V=ifft2(-0.5*(KX.^2+5*KY.^2).*fft2(MinvL0V))-(0.5+mu).*MinvL0V+U.*MinvL0U; 
    MinvL1HermitMinvL0U=ifft2(fft2(L1HermitMinvL0U)./(mu+KX.^2+KY.^2));  
    MinvL1HermitMinvL0V=ifft2(fft2(L1HermitMinvL0V)./(mu+0.5+0.5*KX.^2+2.5*KY.^2));   
        
    if nn == 1
        U=U-MinvL1HermitMinvL0U*DT;
        V=V-MinvL1HermitMinvL0V*DT;
    else
        L1G1=ifft2(-(KX.^2+KY.^2).*fft2(G1))+(V-mu).*G1 + U.*G2; 
        L1G2=ifft2(-0.5*(KX.^2+5*KY.^2).*fft2(G2))-(0.5+mu).*G2 + 0.5*U.*G1; 
        MinvL1G1=ifft2(fft2(L1G1)./(mu+KX.^2+KY.^2)); 
        MinvL1G2=ifft2(fft2(L1G2)./(mu+0.5+0.5*KX.^2+2.5*KY.^2)); 
        MG1=ifft2(fft2(G1).*(mu+KX.^2+KY.^2));  
        MG2=ifft2(fft2(G2).*(mu+0.5+0.5*KX.^2+2.5*KY.^2));     
        alpha=1/sum(sum(MG1.*G1+MG2.*G2))-1/(DT*sum(sum(L1G1.*MinvL1G1+L1G2.*MinvL1G2))); 
        innerproduct=sum(sum(G1.*L1HermitMinvL0U+G2.*L1HermitMinvL0V)); 
        U=U-(MinvL1HermitMinvL0U-alpha*innerproduct*G1)*DT; 
        V=V-(MinvL1HermitMinvL0V-alpha*innerproduct*G2)*DT; 
    end
    
    G1=U-Uold; 
    G2=V-Vold; 
    UVerror(nn)=sqrt(sum(sum(abs(U-Uold).^2+abs(V-Vold).^2)*dx*dy)); UVerror(nn)
    if UVerror(nn) < error_tolerance 
        break
    end
end
   
figure(1)                                   % plotting results
mesh(x, y, U)
figure(2)
mesh(x, y, V)
figure(3)            
semilogy(1:length(UVerror), UVerror, 'linewidth', 2)

