% p6.m: AITEM for computing lattice solitons in the 2D NLS equation
% u_{xx}+u_{yy}-6(sin(x)^2+sin(y)^2)u+u^3=mu*u

 Lx=10*pi; Ly=10*pi; Nx=128; Ny=128; maxerror=1e-10; nmax=500;
 dx=Lx/Nx; x=-Lx/2:dx:Lx/2-dx; kx=[0:Nx/2-1  -Nx/2:-1]*2*pi/Lx;
 dy=Ly/Ny; y=-Ly/2:dy:Ly/2-dy; ky=[0:Ny/2-1  -Ny/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);
 P=1.9092; c=2; DT=0.9;
 u=sech(2*sqrt(X.^2+Y.^2));                 % initial condition
 u=u*sqrt(P/(sum(sum(abs(u.*u)))*dx*dy));   % power normalization
 for nn=1:nmax                              % iteration starts
    L00u=ifft2(-K2.*fft2(u))+(u.*u-V).*u;
    Minvu=ifft2(fft2(u)./(c+K2));
    mu=sum(sum(Minvu.*L00u))/sum(sum(Minvu.*u));
    L0u=L00u-mu*u;
    uerror(nn)=max(max(abs(L0u)));
    if uerror(nn) < maxerror
        break
    end
    u=u+ifft2(fft2(L0u)./(c+K2))*DT;
    u=u*sqrt(P/(sum(sum(abs(u.*u)))*dx*dy));
 end                                        % iteration ends
 mu
 subplot(221); mm=35:94; mesh(x(mm), y(mm), u(mm,mm));
 axis([x(mm(1)) x(mm(end)) y(mm(1)) y(mm(end)) 0 1.3]);
 subplot(222); semilogy(0:nn-1, uerror, 'linewidth', 2);
 axis([0 170 1e-11 10])