%%% On résout l'équation de Fisher-KPP suivante:
%%% u_t = d*u_xx + alpha*u*(1-u/K)

%%% Parameters %%%  
K = 1;
d = 1;
alpha = 1;

%%% KPP nonlinerity 
f = @(u) alpha*u.*(1-u/K);

%%% Discretization of time and space
L = 300;
thresh = 10^(-4);
dt = 0.01; 
dx = 0.03;
N = floor(L/dx);
x = (0:dx:(N-1)*dx)';
Int1 = 0.5;
Int2 = 0.4;
n1 = floor(N*Int1);
n2 = floor(N*Int2);

%%% prepare the operator matrices
E = ones(N,1);
Lap = spdiags([E -2*E E],-1:1,N,N);
Lap(1,1) = -1;
Lap(end,end) = -1;
Lap = Lap/dx^2;
[L_U, U_U] = lu(speye(N)-d*dt*Lap);

%%% default initial conditions
m = floor(0.1*N);
u0 = [ones(m,1);zeros(N-m,1)];
u = u0;

%%% semi-implicit Euler method %%%
nSaveStep=floor(1/dt);
nstep = 1;
figure(1)
while (abs(u(n1))<thresh)
    nstep = nstep+1;
    u = U_U\(L_U\(u + dt*f(u)));
    if (abs(u(n2))<thresh)
        t_p=nstep*dt;
    end
    if (abs(u(n1))<thresh)
        t_m=nstep*dt;
    end
    if mod(nstep,nSaveStep)==0
        plot(x,u,'k','linewidth',2.5) 
        axis([0 L min(u)-0.1 max(abs(u))+0.1]) 
        title(strcat('Time=' ,num2str(nstep*dt)),'fontsize',18)
        drawnow;
   end
end


speed = [N*(Int1-Int2)*dx/(t_m-t_p) 2*sqrt(d*alpha)]
