
%Note: In this case, the solution U for Eq. (134) is real, thus
%   L1= dxx+dyy+V-3U^2+mu,  where V=-V0(sin(x)^2+sin(y)^2). 
%The Hermitian of L1 is itself. 


Lx=10*pi; Ly=10*pi; N=128;                            % mesh 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);

V=-6*(sin(X).^2+sin(Y).^2); P=2.4936; c=1.8; DT=0.6;
U=1.15*sech(X.^2+Y.^2).*cos(X).*cos(Y);               % initial conditions
U=U/sqrt(sum(sum(abs(U).^2))*dx*dy/P);

for nn=1:max_iteration                                % PCSOM iterations start
    Uold=U; 
    L00U=ifft2(-(KX.^2+KY.^2).*fft2(U))+V.*U-U.^3; 
    MinvL00U=ifft2(fft2(L00U)./(c+KX.^2+KY.^2));
    MinvU=ifft2(fft2(U)./(c+KX.^2+KY.^2));
    mu=-sum(sum(U.*MinvL00U))/sum(sum(U.*MinvU));             
    L0U=L00U+mu*U; 
    MinvL0U=ifft2(fft2(L0U)./(KX.^2+KY.^2+c)); 
    L1HermitMinvL0U=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0U))+(V-3*U.^2+mu).*MinvL0U; 
    MinvL1HermitMinvL0U=ifft2(fft2(L1HermitMinvL0U)./(KX.^2+KY.^2+c));
    gamma=sum(sum(U.*MinvL1HermitMinvL0U))/sum(sum(U.*MinvU)); 
    U=U-(MinvL1HermitMinvL0U-gamma*MinvU)*DT; 
    U=U/sqrt(sum(sum(abs(U).^2))*dx*dy/P); 
    Uerror(nn)=sqrt(sum(sum(abs(U-Uold).^2))*dx*dy); Uerror(nn)            
    if Uerror(nn) < error_tolerance 
        break        
    end
end
         
figure(1)                                             % plotting results
mesh(x, y, U)
figure(2)            
semilogy(1:length(Uerror), Uerror, 'linewidth', 2)

