% This code produces Figure S4 in Shaw et al (2023) Ecol Lett
% Pushed to the edge: Spatial Sorting can slow down invasions.

% model written by Frithjof Lutscher

% It plots a population of ten types as it advances in space.
% It computes the speed from the simulation and compares it with
% the analytical expression of the speed in equation (S10)

clear

% parameter definitions
a = 0.02;
mud = 0.05;
p = 0.1:0.1:1;
tau = length(p);
v = 0.25;

% set-up of the spatial grid
l = 50; np = 2^18;  dx = l/np;  x = linspace(-l/2,l/2-dx,np);    
xx = linspace(-l, l-dx, 2*np);   % double domain for padding
PAD = (abs(xx)<=l/2); % padding outside [-l/2, l/2]

% Dispersal kernel with parameter scaled to 1
Lap = exp(-sqrt(2/v)*abs(xx))/sqrt(2*v);
FLap = fft(Lap);

% initial conditions
for i = 1:tau
n(i,:) = 1/tau*(abs(xx)<1);
end

% The iteration loop for 5 time steps
for kk=1:5
  Ntot = sum(n);  % total population
  % the growth step
  CUT = (Ntot>a);
  extent(kk) = max(find(Ntot>a))*dx-l;
  for i=1:tau
    fn(i,:) = 1/tau*CUT;   
  end
  % the dispersal step
  for i = 1:tau
  nhat(i,:) = p(i)*(1-mud)*dx*real( fftshift( ifft( fft(fn(i,:)).*FLap ) ) );
  nhat(i,:) = nhat(i,:) + (1-p(i))*fn(i,:);
  end
 
% plotting the figure 
  figure(1)
  for i = 1:tau
    plot(xx,n(i,:),'linewidth',2)
    hold on 
    plot(xx,fn(i,:))
    plot(xx,nhat(i,:))
  end
  hold off
  axis([0 6 0 0.15])
  pause(0.01)
  n = nhat;
  clear fn nhat

end

% Now calculate the approximate speed from the simulation by taking the difference
% in extent for the verious times
c_sim = diff(extent)
c_sim(end)

% the analytic formula in equation (S10)
b = sqrt(v/2);
c = b*log((1-mud)*sum(p)/(2*a*tau))