% p8.m: the Newton-CG method 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=20*pi; Ly=20*pi; N=256; maxerror=1e-10; errorCG=1e-2;
 x=-Lx/2:Lx/N:Lx/2-Lx/N; kx=[0:N/2-1  -N/2:-1]*2*pi/Lx;
 y=-Ly/2:Ly/N:Ly/2-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); c=3; mu=-4.56;
 u=0.56*sech(0.5*sqrt(X.^2+Y.^2)).*cos(X).*cos(Y);
 nnt=0;                  % nnt: number of Newton iterations
 ncg=0;                  % ncg: total number of CG iterations
 while 1                 % Newton iterations begin
    nnt=nnt+1;
    L0u=ifft2(-K2.*fft2(u))-(V+u.*u+mu).*u;
    uerror(nnt)=max(max(abs(L0u))); uerror(nnt)
    numcg(nnt)=ncg;
    if uerror(nnt) < maxerror
        break
    end
    du=0*u;                     % CG iterations begin
    R=-L0u; MinvR=ifft2(fft2(R)./(c+K2));
    R2=sum(sum(R.*MinvR)); R20=R2;
    D=MinvR;
    while (R2 > R20*errorCG^2)
        L1D=ifft2(-K2.*fft2(D))-(V+3*u.*u+mu).*D;
        a=R2/sum(sum(D.*L1D));
        du=du+a*D;
        R=R-a*L1D; MinvR=ifft2(fft2(R)./(c+K2));
        R2old=R2;
        R2=sum(sum(R.*MinvR));
        b=R2/R2old;
        D=MinvR+b*D;
        ncg=ncg+1;
    end                         % CG iterations end
    u=u+du;
 end                     % Newton iterations end
 subplot(221);           % solution plotting
 nr=70:1:188; mesh(x(nr),y(nr),u(nr, nr));
 axis([x(70) x(188) x(70) x(188) -0.3 0.5]); view([-20 20]);
 xlabel('x'); ylabel('y'); zlabel('u', 'rotation',0);
 subplot(222);
 semilogy(numcg, uerror, numcg, uerror, 'o');
 axis([0 300 1e-12 1]);
 xlabel('number of iterations'); ylabel('error')