%%   Simulate Fisher's equation in 1D

%%%%%%%%%%%%%%%%% Numerical/problem parameters
clear;  %  gets rid of old data
L = 15;  % domain size
dt = 0.2; % time step size
N = 300; % number of spatial mesh points
dx = L/N;
x = [dx: dx: L-dx]';  %setup grid
D = .01;     %Diffusion coefficient

%%%%%%%%%%%%%%%%% Nonlinearity (comment/uncomment this)

f = inline('u.*(1-u)','u');   %logistic growth

%f = inline('u.*(u-.4).*(1-u)','u');   %logistic growth with harvesting

%%%%%%%%%%%%%%%%%%%%%   Initial data (comment/uncomment this)

% "seed" with sharp edge
seedsize = 1.0;
u = (x>3).*(x<3+seedsize);

% initial seed with less spatial decay (leads to faster propagation?)
decay = .05;  %decay constant 
%u = exp( -abs(x-3.5)/decay);

%%%%%%%%%%%%%%%%%%%%% Other stuff

%diffusion matrix
s = D*dt/dx^2; 
A = (1+2*s)*diag(ones(N-1,1),0) - s*diag(ones(N-2,1),1) - s*diag(ones(N-2,1),-1);
A(1,1) = A(1,1) - s;
A(end,end) = A(end,end)-s;   % neumann boundary conditions

%%%%%%%%%%%%%%%%%%%%% Time stepping
time=0;  
while(1),
   plot(x,u);
   text(L/2,1,['time=' num2str(time)]);
   axis([0 L 0 2])
   xlabel('x');  ylabel(' u(x,t)');
   pause(.05);
   unext = A\u;        % Solve the equation A*unext =u;
   u = unext + dt*f(u);    % now u is equal to the updated solution
   time = time + dt;
end;
