%Note 1:  Let U=u+iv. The equations for (u, v) are Eqs. (135)-(136) in the paper. %       Operator L1 = [dxx+dyy+V+3u^2+v^2+mu           2uv%                             2uv              dxx+dyy+V+3v^2+u^2+mu],%       where V=-V0(sin(x)^2+sin(y)^2). %   Hermitina of L1 is L1 itself. %Note 2:  In the computer execution, equations for u and v are recombined %into one for U=u+iv for numerical efficiency. In doing so, for F=f+ig, % L1[f, g]' is recombined into  Fxx+Fyy+(V+|U|^2+mu)F+2U*real(conj(U)*F)%Note 3: For complex numbers A and B, real(A)*real(B)+imag(A)*imag(B)=real[conj(A)B] %This fact is used in the calculations of mu and gamma in lines 30 and 35Lx=12*pi; Ly=12*pi; N=256;                              % mesh parametersmax_iteration=1e5; 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=14.6004; c=3.7; DT=0.8;   % other parameters and initial conditionsU=1.7*(exp(-X.^2-Y.^2)+i*exp(-(X-pi).^2-Y.^2)-exp(-(X-pi).^2-(Y-pi).^2)-i*exp(-X.^2-(Y-pi).^2)); 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+abs(U).^2).*U;     MinvL00U=ifft2(fft2(L00U)./(c+KX.^2+KY.^2));    MinvU=ifft2(fft2(U)./(c+KX.^2+KY.^2));    mu=-sum(sum(real(conj(U).*MinvL00U)))/sum(sum(conj(U).*MinvU));        L0U=L00U+mu*U;     MinvL0U=ifft2(fft2(L0U)./(KX.^2+KY.^2+c));     L1HermitMinvL0U=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0U))+(V+abs(U).^2+mu).*MinvL0U+2*U.*real(conj(U).*MinvL0U);     MinvL1HermitMinvL0U=ifft2(fft2(L1HermitMinvL0U)./(KX.^2+KY.^2+c));    gamma=sum(sum(real(conj(U).*MinvL1HermitMinvL0U)))/sum(sum(conj(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    endendfigure(1)                                                % plotting resultsimagesc(x, y, abs(U))axis([-2.5*pi 3.5*pi -2.5*pi 3.5*pi])figure(2)imagesc(x, y, angle(U))axis([-2.5*pi 3.5*pi -2.5*pi 3.5*pi]) figure(3)            semilogy(1:length(Uerror), Uerror, 'linewidth', 2)