% [u,x,t,kk,amp,M,H] = split_nonlinear_heat(t_0,t_f,dt,NPrint,N)
%
% this solves the nonlinear heat equation 
%
%          u_t = u_xx + u^2
%
% with periodic boundary conditions using spectral derivatives in
% space and a splitting method for the time stepping.  t_0 is the 
% initial time, t_f is the
% final time, dt is the timestep size, N is the number of mesh-points, and 
% NPrint is the number of printouts.  It also returns the power spectrum of the
% solution, where kk is the collection of wave numbers and amp = log(f_k^2)
% is the magnitude of the kth Fourier mode.
%
% N SHOULD BE A POWER OF 2!  
%
% (It will work no matter what, but will go much faster if N is a power
% of 2.)

function [u,x,t,kk,amp] = split_nonlinear_heat(t_0,t_f,dt,NPrint,N)

% define the mesh in space
dx = 2*pi/N;
x = 0:dx:2*pi-dx;
x = x';
% the wave-numbers...
kk = (-N/2+1):1:(N/2-1);

% total number of steps to be taken:
Ntot = (t_f-t_0)/dt;
% test to see if Ntot is an integer, to see if your choice of dt
% is compatible with your choice of t_0 and t_f
if abs(Ntot - round(Ntot)) > 100000*eps
  disp(' ')
  disp('dt is not consistent with your choice of t_0 and t_f.  Choose another dt!')
  disp(' ')
  break
end
Ntot = round(Ntot);
% number of steps to be taken before a (screen) printout
Nrun = Ntot/NPrint;
% test to see if Nrun is an integer to see if your choice of dt is
% consistent with NPrint
if abs(Nrun - round(Nrun)) > 1000*eps
  disp(' ')
  disp('NPrint is not consistent with dt.  Choose another NPrint!')
  disp(' ')
  break
end
Nrun = round(Nrun);
% define the vector t that is to be printed out...
t = t_0:Nrun*dt:t_f;

% choose initial data
u_now = cos(3*x);

% We're computing solutions of
%
%    u_t = u_xx + |u|^2
%
% the first step is to take one time-step of size dt
% to compute the solution of u_t = u_xx
% To do this, take the FFT of u, resulting in 
%
%        u_k exp(i k x)
%
% We know how to solve u_t = u_xx at the spectral level,
% resulting in 
%
%       u_k exp(- |k|^2 dt) exp(i k x)
%
% Taking the IFFT, we have a function v
%
%   v = IFFT( exp(- |k|^2 dt) FFT(u) )
% 
% We now use v as initial data to take one timestep of size
% dt to sove the problem u_t = u^2
% this is just an ODE and at each meshhpoint it has the solution
%
%    u(x_j)/(1 - dt u(x_j))
%
% We take this function as our solution u at time t + dt to
% the NLS.

% save first printout:
u(:,1) = u_now;
% save the power spectrum 
v = fft(u_now);
amp(N/2:N-1,1) = log10( abs(v(1:N/2)) + 1.e-16);
amp(1:N/2-1,1) = log10( abs(v(N/2+2:N)) + 1.e-16);
for jj = 1:NPrint
  for j=1:Nrun

    % advance with with one step of size dt
    u_coarse = tstep(u_now,dt,N);
    % advance with with two steps of size dt/2
    u_fine = tstep(u_now,dt/2,N);
    u_fine = tstep(u_fine,dt/2,N);
    % use Richardson extrapolation to combine u_coarse and
    % u_fine to get a solution with local truncation
    % error of O(dt^3)...
    u_now = 2*u_fine - u_coarse;

  end
  % save printout
  u(:,jj+1) = u_now;
  v = fft(u_now);
  amp(N/2:N-1,jj+1) = log10( abs(v(1:N/2)) + 1.e-16);
  amp(1:N/2-1,jj+1) = log10( abs(v(N/2+2:N)) + 1.e-16);
end
x(N+1) = 2*pi;
u(N+1,:) = u(1,:);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function u_now = tstep(u_now,dt,N);

% first take one time-step of the equation u_t = u_xx
v = fft(u_now);
ii = 1:N/2;
k = (ii-1)';
v(ii) = exp(-k.^2*dt).*v(ii);    
ii = N/2+2:N;
k = (ii-N-1)';
v(ii) = exp(-k.^2*dt).*v(ii);    
v = real(ifft(v));

% now take one time-step of u_t = u^2
% using v as our initial data
u_now = v./(1-dt*v);