% This code uses Matlab's syntax.
 % This file is best viewed with Notepad or Wordpad, 
 %  but not with Notepad++ (some of the lines will appear as shifted there).

 % This is a logical flow-chart for implementing the averaging equation, Eq. (28),
 % of the paper "Long-time simulations of nonlinear Schr¨odinger-type equations 
 % using step size exceeding threshold of numerical instability," by T.I. Lakoba, 
 % J. Sci. Comput.  DOI 10.1007/s10915-016-0346-y
 %
 % Recall from the line below Eq. (25) that 
 %       \hat{v}(n) = \hat{u}(n) * exp(i*k^2*dt*n).              		(1)
 %  Here \hat{u}(n) is the Fourier transform of the solution u at time step n, 
 %  and dt is the time step. 
 %  For the purpose of coding, it is convenient to write (28) as:
 % 		 vhat_new = c(1)*vhat(1) + c(2)*vhat(2) + ... c(7)*vhat(7),	 (2)
 %  where the averaging coefficients  c(1) -- c(7) are defined below, and 
 %  vhat(1) stands for \hat{v}(n-3), vhat(2) stands for \hat{v}(n-2), etc..
 %  Note:   The variables \hat{v}(n) depend not only on time (=dt*n) but also on wavenumber, k.
 %          However, since this is text is just the logical flow-chart rather than an actual code,
 %          we will omit the k-dependence (and also the x-dependence later). 
 %  Finally, note that equation (2) above can be rewritten in terms of variables \hat{u}:
 %       uhat_new = c(1)*uhat(1)*exp(-3*i*k^2*dt) + c(2)*uhat(2)*exp(-2*i*k^2*dt) + ... ,
 %                                                                                (3)
 %  where we have used equation (1) relating vhat and uhat. 
 

 % In the beginning of the code we specify how often we apply stabilization step:
        tstabilize = 5; 	% time interval over which we stabilize the solution
	nn_tstabilize = round(tstabilize/dt);
				% the corresponding number of time steps
	  
 % Coefficients needed for averaging, as per Eq. (28):
        c(1) = 1/64;
	c(2) = -3/32;
	c(3) = 15/64;
	c(4) = 11/16;
	c(5) = c(3);
	c(6) = c(2);
	c(7) = c(1);
	 
 % We also define exponential factors to relate \hat{v} with \hat{u}, as stated above:
        ef(1) = exp(-3*i*k^2*dt); 
	ef(2) = exp(-2*i*k^2*dt); 
	ef(3) = exp(-1*i*k^2*dt); 
	ef(4) = 1; 
	ef(5) = exp(1*i*k^2*dt); 
	ef(6) = exp(2*i*k^2*dt); 
	ef(7) = exp(3*i*k^2*dt); 
 % As with vhat(n), these exponential factors depend on k. However, we do not indicate
 %  such a dependence explicitly in this logical flow-chart. 
 
	  
  % -------------------------
  % Calculations begin here:
     
for n = 1 : nmax 
	 
	u = ...     % this is solution u(n) computed by Leap-frog
	u1 = ...    % this is solution u(n+1) computed by Leap-frog
		
	if  n > 4   % we simply want the stabilization step being applied at very first steps
			  
		% We collect 8 = (7+1) values of u to use in the averaging formula (28).
		%  The reason why we need 1 more value than in (28) will become apparent soon:
		%  we need to obtain a new value not only for u(n) but also for u(n+1).
		for j = 1 : 8
			if  mod(n + 3, nn_tstabilize) == j-1 
				ustab(j) = u;     % this is the collected value of the solution;
			      	           	  % recall that u also depends on x, but we do not indicate
						  %  such a dependence in this flow-chart
			end
		end
				  
     		if  mod(n, nn_tstabilize) == 4
			% In order to avoid the accummulation of the parasitic error, we
			%  constract both u(n) and u(n+1) from scratch.
			u = 0;
			u1 = 0;
				  
			% "Build" the new solutions u(n) and u(n+1) by adding one term at a time
			%  as per Eq. (28).
			for j = 1 : 7
				u = u + c(j)*ifft( ef(j).* fft(ustab(j)) );
				u1 = u1 + c(j)*ifft( ef(j).*  fft(ustab(j+1)) );
			end
			% Note 1:
			%  For those unfamiliar with Matlab's syntax, fft and ifft stand for the
               		%  Fast Fourier Transform and its inverse.
                	% Note 2:
                	%   The occurrence of the exponential factors is clarified by equation (3) above.				
           	 end
		 % Strictly speaking, above we have reset the u,u1-values 3 time steps back and so 
		 % need to also reset the time (i.e. the step counter n) 3 steps back.
         	 % However, not doing so will introduce a negligible error (proportional to 
		 %  3/nn_tstabilize << 1) in the total simulated time, and so we won't do it
		 %  in order not to make the logic of the code more complex.
			
	end		% this ends the loop verifying if n > 4
				  
end             	% this ends the loop of computational steps