
%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);  

P=47.3744;   % this power is the integral of u^2+v^2
DT=0.63; 
U=sech(0.35*sqrt(X.^2+Y.^2));               % initial conditions
V=0.3*sech(0.5*sqrt(X.^2+Y.^2)); 
U=U/sqrt(sum(sum(U.^2+V.^2))*dx*dy/P);
V=V/sqrt(sum(sum(U.^2+V.^2))*dx*dy/P);

for nn=1:max_iteration                      % PCSOM iterations start
    Uold=U; 
    Vold=V; 
    L00U=ifft2(-(KX.^2+KY.^2).*fft2(U))+U.*V; 
    L00V=0.5*(ifft2(-(KX.^2+5*KY.^2).*fft2(V))+0.5*U.^2-V); 
    if nn == 1
        mu=sum(sum(U.*L00U+V.*L00V))/sum(sum(U.^2+V.^2));
    end
    MinvL00U=ifft2(fft2(L00U)./(mu+KX.^2+KY.^2));  
    MinvL00V=ifft2(fft2(L00V)./(mu+0.5+0.5*KX.^2+2.5*KY.^2));  
    MinvU=ifft2(fft2(U)./(mu+KX.^2+KY.^2));  
    MinvV=ifft2(fft2(V)./(mu+0.5+0.5*KX.^2+2.5*KY.^2));  
    mu=sum(sum(U.*MinvL00U+V.*MinvL00V))/sum(sum(U.*MinvU+V.*MinvV));
    L0U=L00U-mu*U; 
    L0V=L00V-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));   
    gamma=sum(sum(U.*MinvL1HermitMinvL0U+V.*MinvL1HermitMinvL0V))/sum(sum(U.*MinvU+V.*MinvV)); 
	U=U-(MinvL1HermitMinvL0U-gamma*MinvU)*DT;
    V=V-(MinvL1HermitMinvL0V-gamma*MinvV)*DT;
    U=U/sqrt(sum(sum(U.^2+V.^2))*dx*dy/P);
    V=V/sqrt(sum(sum(U.^2+V.^2))*dx*dy/P);
 
    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)














