% fisher.m 1D model of Fisher's equation
% variable is u = N/N(max)
% only scaling variable is D
npoints = 80; % Number of points along x-axis 
start = 20; % Start with a spike at this location
u = zeros(1,npoints);
u = -ones(1,npoints);
u(1:start) = 5;
tmax = 200; 
dt = 0.2;
dx = 0.5;
xmax = npoints*dx;
x = dx:dx:xmax;

D=0.2; % Diffusion constant
r = 0.5; % base rate of growth
epsilon = 0.1;

figure(1)
uplot = -1:0.01:4;
dudt = -r*uplot.*(1-uplot)-epsilon*uplot.*uplot.*uplot;

plot(uplot,dudt,'r');
hold on
plot(uplot,zeros(size(uplot)));
xlabel('u')
ylabel('dudt')

f = D*dt/(dx*dx)
%% Now integrate through time
for t = 0:dt:tmax; 
    % Next the diffusion term, calculating second derivative 
    u(2:end-1) = (1-f)*u(2:end-1)+0.5*f*(u(3:end) + u(1:end-2));
    
    % Care for boundary effects
    u(1) = (1-f)*u(1) + 0.5*f*u(2);
    u(end) = (1-f)*u(end)+0.5*f*u(end-1);

    % Next the population growth term
    u = u + (-r*u.*(1-u)-epsilon*u.*u.*u)*dt;

    figure(2);
    plot(x,u);
    axis([0 xmax -1 4])
    xlabel('position')
    ylabel('u')
end
