% p7.m: the MSOM for computing gap solitons in the 2D NLS equation
% with a periodic potential and self-defocusing nonlinearity:
% u_{xx}+u_{yy}-6(sin(x)^2+sin(y)^2)u-u^3=mu*u

 Lx=10*pi; Ly=10*pi; N=128; nmax=1000; maxerror=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); K2=KX.^2+KY.^2;
 V=6*(sin(X).^2+sin(Y).^2);
 mu=-5; u=1.15*sech(sqrt(X.^2+Y.^2)).*cos(X).*cos(Y);
 c=2.9; DT=1.7;
 for nn=1:nmax                         % MSOM iterations start
    uold=u;
    F0=-V-u.*u-mu; F1=-V-3*u.*u-mu;
    L0u=ifft2(-K2.*fft2(u))+F0.*u;
    uerror(nn)=max(max(abs(L0u))); uerror(nn)
    if uerror(nn) < maxerror
        break
    end
    MinvL0u=ifft2(fft2(L0u)./(c+K2));
    L1MinvL0u=ifft2(-K2.*fft2(MinvL0u))+F1.*MinvL0u;
    MinvL1MinvL0u=ifft2(fft2(L1MinvL0u)./(c+K2));
    if nn == 1
        u=u-MinvL1MinvL0u*DT;
    else
        L1G=ifft2(-K2.*fft2(G))+F1.*G;
        MinvL1G=ifft2(fft2(L1G)./(c+K2));
        MG=ifft2(fft2(G).*(c+K2));
        alpha=1/sum(sum(MG.*G))-1/(DT*sum(sum(L1G.*MinvL1G)));
        u=u-(MinvL1MinvL0u-alpha*sum(sum(G.*L1MinvL0u))*G)*DT;
    end
    G=u-uold;
 end                                   % MSOM iterations end
 subplot(221); mm=35:94; mesh(x(mm), y(mm), u(mm,mm));
 axis([x(35) x(94) y(35) y(94) 0 1.3])
 subplot(222); semilogy(0:nn-1, uerror, 'linewidth', 2);
 axis([0 350 1e-11 10])