
Lx=20; N=128;                                      % grid parameters
x=-Lx/2:Lx/N:Lx/2-Lx/N;          
kx=[0:N/2-1  -N/2:-1]*2*pi/Lx; 
max_iteration=1e4; error_tolerance=1e-10; 
mu=1; B=1/sqrt(1+16*mu/3);                         % soliton solution
u=sqrt(4*B*mu./(cosh(2*sqrt(mu)*x)-B)); 
c=2.4; DT=2.3;                                     % scheme parameters
F1=-mu-u.^2+u.^4; 
F2=-mu-3*u.^2+5*u.^4; 
V=-0.3*(1-2*x.^2).*exp(-x.^2);                     % initial conditions
W=i*exp(-x.^2);      

for nn=1:max_iteration                             % OOM iterations begin
    % first RK4 step
    LV=-i*(ifft(-kx.^2.*fft(W))+F1.*W); 
    LW=-i*(ifft(-kx.^2.*fft(V))+F2.*V);
    MinvV=ifft(fft(V)./(c+kx.^2));
    MinvW=ifft(fft(W)./(c+kx.^2));
    lambda=sum(conj(MinvV).*LV+conj(MinvW).*LW)   ...
          /sum(conj(MinvV).*V+conj(MinvW).*W);
    Lambda(nn)=lambda;
    K1V=ifft(fft(LV-lambda*V)./(c+kx.^2));  
    K1W=ifft(fft(LW-lambda*W)./(c+kx.^2));  
    TV=V+0.5*K1V*DT; 
    TW=W+0.5*K1W*DT; 
    % second RK4 step
    LTV=-i*(ifft(-kx.^2.*fft(TW))+F1.*TW); 
    LTW=-i*(ifft(-kx.^2.*fft(TV))+F2.*TV);
    MinvTV=ifft(fft(TV)./(c+kx.^2));
    MinvTW=ifft(fft(TW)./(c+kx.^2));
    Tlambda=sum(conj(MinvTV).*LTV+conj(MinvTW).*LTW)  ...
           /sum(conj(MinvTV).*TV+conj(MinvTW).*TW);
    K2V=ifft(fft(LTV-Tlambda*TV)./(c+kx.^2));  
    K2W=ifft(fft(LTW-Tlambda*TW)./(c+kx.^2));  
    TV=V+0.5*K2V*DT; 
    TW=W+0.5*K2W*DT; 
    % third RK4 step
    LTV=-i*(ifft(-kx.^2.*fft(TW))+F1.*TW); 
    LTW=-i*(ifft(-kx.^2.*fft(TV))+F2.*TV);
    MinvTV=ifft(fft(TV)./(c+kx.^2));
    MinvTW=ifft(fft(TW)./(c+kx.^2));
    Tlambda=sum(conj(MinvTV).*LTV+conj(MinvTW).*LTW)  ...
           /sum(conj(MinvTV).*TV+conj(MinvTW).*TW);
    K3V=ifft(fft(LTV-Tlambda*TV)./(c+kx.^2));  
    K3W=ifft(fft(LTW-Tlambda*TW)./(c+kx.^2));  
    TV=V+K3V*DT; 
    TW=W+K3W*DT; 
    % fourth RK4 step
    LTV=-i*(ifft(-kx.^2.*fft(TW))+F1.*TW); 
    LTW=-i*(ifft(-kx.^2.*fft(TV))+F2.*TV);
    MinvTV=ifft(fft(TV)./(c+kx.^2));
    MinvTW=ifft(fft(TW)./(c+kx.^2));
    Tlambda=sum(conj(MinvTV).*LTV+conj(MinvTW).*LTW)  ...
           /sum(conj(MinvTV).*TV+conj(MinvTW).*TW);
    K4V=ifft(fft(LTV-Tlambda*TV)./(c+kx.^2));  
    K4W=ifft(fft(LTW-Tlambda*TW)./(c+kx.^2));  
    V=V+(K1V+2*K2V+2*K3V+K4V)*DT/6; 
    W=W+(K1W+2*K2W+2*K3W+K4W)*DT/6; 
    
    if nn == 1
        lambda_error=1; 
        else
        lambda_error=abs(Lambda(nn)-Lambda(nn-1))
    end
    if lambda_error < error_tolerance
        break
    end
end
lambda
semilogy(abs(Lambda(2:nn)-Lambda(1:nn-1)))         % plotting error diagram

