% p5.m: Petviashvili method for computing the ground state in the
% 2D NLS equation: u_{xx}+u_{yy}+u^3=u.

 Lx=30; Ly=30; Nx=128; Ny=128; maxerror=1e-10; nmax=300;
 x=-Lx/2:Lx/Nx:(Lx/2-Lx/Nx); y=-Ly/2:Ly/Ny:(Ly/2-Ly/Ny);
 kx=[0:Nx/2-1 -Nx/2:-1]*2*pi/Lx; 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;
 p=3; gamma=p/(p-1);
 mu=1; u=2.2*exp(-(X.^2+Y.^2));
 for nn=1:nmax
    Mu=ifft2((mu+K2).*fft2(u));
    u3=u.*u.*u;
    L0u=u3-Mu;
    uerror(nn)=max(max(abs(L0u)));
    if uerror(nn) < maxerror
        break
    end
    S=sum(sum(u.*Mu))/sum(sum(u.*u3));
    u=S^gamma*ifft2(fft2(u3)./(mu+K2));
 end
 subplot(221)
 mm=36:94; mesh(x(mm), y(mm), u(mm,mm));
 axis([x(mm(1)) x(mm(end)) y(mm(1)) y(mm(end)) 0 2]);
 xlabel('x'); ylabel('y'); zlabel('|u|'); title('(a)');
 subplot(222)
 semilogy(0:nn-1, uerror, 'linewidth', 2);
 axis([0 50 1e-11 10]);
 xlabel('number of iterations'); ylabel('error'); title('(b)')