Lx=10*pi;  Ly=10*pi;  N=128;  c=0.7;  DT=1;
max_iteration=10000; error_tolerance=1e-10;

dx=Lx/N; x=-Lx/2:dx:Lx/2-dx; kx=[0:N/2-1 -N/2:-1]*2*pi/Lx;
dy=Ly/N; y=-Ly/2:dy:Ly/2-dy; 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=3*(cos(X).^2+cos(Y).^2); 

P=3;
U=exp(-(X.^2+Y.^2));  
U=U/sqrt(sum(sum(abs(U).^2))*dx*dy/P);
for nn=1:max_iteration
    Uold=U;
    L00U=ifft2(-K2.*fft2(U))+V.*U + U.^3;
    MinvU=ifft2(fft2(U)./(K2+c));
    mu=sum(sum(L00U.*MinvU))/sum(sum(U.*MinvU));
    U=U+ifft2(fft2(L00U-mu*U)./(K2+c))*DT;
    U=U/sqrt(sum(sum(abs(U).^2))*dx*dy/P);
    Uerror(nn)=sqrt(sum(sum(abs(U-Uold).^2))*dx*dy); 
    Uerror(nn)
    if Uerror(nn) < error_tolerance
        break
    end
end
mesh(x, y, U)