 Lx=20*pi; Ly=20*pi; N=256; errormax=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);

 ncg=0; nnt=0;
 while 1
    nnt=nnt+1;
    L0U=ifft2(-K2.*fft2(U))+(V-U.*U+mu).*U;
    errorU(nnt)=max(max(abs(L0U)));
    iterations(nnt)=ncg;
    errorU(nnt)
    if errorU(nnt) < errormax
        break
    end
    DU=0*X;
    R=-L0U; MinvR=ifft2(fft2(R)./(K2+c));
    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)./(K2+c));
        R2old=R2;
        R2=sum(sum(R.*MinvR));
        b=R2/R2old;
        D=MinvR+b*D;
        ncg=ncg+1;
    end
    U=U+DU;
 end
figure(1); nr=60:1:198; mesh(x(nr),y(nr),U(nr, nr))
axis([x(60) x(198) x(60) x(198) -0.3 0.5]); view([-20 20])
xlabel('x', 'fontsize', 18); ylabel('y', 'fontsize', 18); zlabel('|u|', 'fontsize', 18)
figure(2); semilogy(iterations, errorU, iterations, errorU, 'ro'); 
xlabel('number of CG iterations', 'fontsize', 16); ylabel('error', 'fontsize', 16)


