
%Note: For Eqs. (146)-(147), 
% L1 = [dxx+dyy+T11-mu1        T12
%           T12          dxx+dyy+T22-mu2]
%where T11=(3*U^2+V^2+s*(U^2+V^2)^2)/(1+s*(U^2+V^2))^2; 
%      T22=(3*V^2+U^2+s*(U^2+V^2)^2)/(1+s*(U^2+V^2))^2; 
%      T12=2UV/(1+s(U^2+V^2))^2
%The Hermitian of L1 is itself


Lx=24; Ly=24; N=64;                          % mesh parameters
max_iteration=10000; error_tolerance=1e-10; 
x=-Lx/2:Lx/N:Lx/2-Lx/N; dx=Lx/N; kx=[0:N/2-1  -N/2:-1]*2*pi/Lx;
y=-Ly/2:Ly/N:Ly/2-Ly/N; dy=Ly/N; ky=[0:N/2-1  -N/2:-1]*2*pi/Ly;
[X, Y]=meshgrid(x, y); [KX, KY]=meshgrid(kx, ky); [TH R]=cart2pol(X,Y); 

s=0.5; mu1=1; mu2=0.5; DT=1.9; 
U=3*exp(-0.2*R.^2);                          % initial conditions
V=1.5*R.*exp(-0.2*R.^2).*cos(TH); 

for nn=1:max_iteration                       % SOM iterations start
	Uold=U; 
    Vold=V; 
    L0U=ifft2(-(KX.^2+KY.^2).*fft2(U))+(U.^2+V.^2)./(1+s*(U.^2+V.^2)).*U-mu1*U; 
    L0V=ifft2(-(KX.^2+KY.^2).*fft2(V))+(U.^2+V.^2)./(1+s*(U.^2+V.^2)).*V-mu2*V; 
    MinvL0U=ifft2(fft2(L0U)./(mu1+KX.^2+KY.^2));  
    MinvL0V=ifft2(fft2(L0V)./(mu2+KX.^2+KY.^2));        
    T11=(3*U.^2+V.^2+s*(U.^2+V.^2).^2)./(1+s*(U.^2+V.^2)).^2; 
    T22=(3*V.^2+U.^2+s*(U.^2+V.^2).^2)./(1+s*(U.^2+V.^2)).^2; 
    T12=2*U.*V./(1+s*(U.^2+V.^2)).^2;  
    L1HermitMinvL0U=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0U))+(T11-mu1).*MinvL0U + T12.*MinvL0V; 
    L1HermitMinvL0V=ifft2(-(KX.^2+KY.^2).*fft2(MinvL0V))+(T22-mu2).*MinvL0V + T12.*MinvL0U; 
    MinvL1HermitMinvL0U=ifft2(fft2(L1HermitMinvL0U)./(mu1+KX.^2+KY.^2));  
    MinvL1HermitMinvL0V=ifft2(fft2(L1HermitMinvL0V)./(mu2+KX.^2+KY.^2));   
    U=U-MinvL1HermitMinvL0U*DT;
    V=V-MinvL1HermitMinvL0V*DT;
    G1=U-Uold; 
    G2=V-Vold; 
    UVerror(nn)=sqrt(sum(sum(abs(U-Uold).^2+abs(V-Vold).^2)*dx*dy)); UVerror(nn)
    if UVerror(nn) < error_tolerance 
	    break
    end
end

figure(1)                                      % plotting results
mesh(x, y, U)
figure(2)
mesh(x, y, V)
figure(3)            
semilogy(1:length(UVerror), UVerror, 'linewidth', 2)

