% simulate the dispersal renewal process

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

 

% set up the grid

L = 100; % size of domain
N = 1000; % number of grid points
h = L/N;
x = [0:N]*h;
XY = ones(N+1,1)*x-x'*ones(1,N+1);

bet = 2; %this is the drift velocity  (see exercise 6)
bet = 0;
bsq = sqrt(bet^2+4);

kp = (bsq+bet)/2;
km = (bsq-bet)/2;

% the dispersal kernal
 
K = exp(-(bsq*abs(XY)-bet*XY)/2)/bsq;
A = 0.5;
%initialize u

u =0.5*sech(0.5*(x'-L/2));


figure(2)
plot(x,u,'--','linewidth',2)
hold on
U = u;

% now iterate 
for j = 1:35
F = u ./(A+u );  % here is where we specify F
 
 u = K*F*h;  %This is the iteration
 U = [U,u];  % save the answer
 
end
figure(2)
 plot(x,U) % plot the answer
  
xlabel('x','fontsize',20)
ylabel('u_n(x)','fontsize',20)
hold off
 
 
figure(3)
lvl=(1-A)
contour(U,[lvl lvl],'linewidth',2)

xlabel('n','fontsize',20)
ylabel('x','fontsize',20)

    
