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

figure(1)
figure(2)
clf

%%

%%% Parameters %%%  
K = 1;
d = 1;
alpha = 1;
c = 0.5; %%% speed of climate shift
l = 3; %%% length of favorable region
r = 1; %%% death rate outside favorable region
Xp = -c/2+sqrt(1+(c/2)^2);
Xm = -c/2-sqrt(1+(c/2)^2);


%%% Discretization of time and space
L = 200;
thresh = 10^(-4);
dt = 0.01; 
dx = 0.03;
N = floor(2*L/dx)+1;
x = (-L:dx:L)';
Int1 = 0.9;
Int2 = 0.85;
n1 = floor(N*Int1);
n2 = floor(N*Int2);
nstop = floor(3*N/4);

%%% Heterogeneous nonlinerity 
g = @(u,y) alpha*(1-u/K).*(0<y).*(y<l)-r*(y<0)-r*(y>l);


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

%%% 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);

%% Résolution

%%% semi-implicit Euler method %%%
nSaveStep=floor(1/dt);
nstep = 1;
t = 0;
figure(1)
while (abs(u(nstop))<thresh)
    nstep = nstep+1;
    u = U_U\(L_U\(u + dt*u.*g(u,x-c*t)));
    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
   t = t+dt;
end

clf

%%     

Ind = find(u>thresh);
U = u(Ind);
V = diff(U)/dx;
U = U(2:end);

I1 = find(Xp*U-V>0.01);
N1 = Ind(I1(1))-1;
I2 = find(V-Xm*U<0);
N2 = Ind(I2(1));

figure(1)
axes('linewidth',2,'FontWeight','bold')
hold
plot(x(1:N1)-x(N1),u(1:N1),'k','linewidth',2.5)
plot(x(N1:N2)-x(N1),u(N1:N2),'b','linewidth',2.5)
plot(x(N2:end)-x(N1),u(N2:end),'k','linewidth',2.5)
plot(zeros(length(x),1),x,'k--','linewidth',1)
plot(l*ones(length(x),1),x,'k--','linewidth',1)
axis([-15 15 min(u)-0.05 max(u)+0.05])
xlabel('x','FontSize',16,'FontWeight','bold')
ylabel('u','FontSize',16,'FontWeight','bold')
box on


%% Portrait de phase
v = diff(u)/dx;
u = u(2:end);
figure(2)
axes('linewidth',2,'FontWeight','bold')
hold on
plot(x,Xp*x,'r','linewidth',1.5)
plot(x,Xm*x,'r','linewidth',1.5)
plot(x,zeros(length(x),1),'k','linewidth',1.25)
plot(zeros(length(x),1),x,'k','linewidth',1.25)
plot(u(N1:N2),v(N1:N2),'b','linewidth',2.5)
plot(u(1:N1),v(1:N1),'k','linewidth',3.5)
plot(u(N2:end),v(N2:end),'k','linewidth',3.5)
axis([-0.1 1.1 -0.75 0.75])
xlabel('u','FontSize',16,'FontWeight','bold')
ylabel('v','FontSize',16,'FontWeight','bold')
box on

% %% Wave speed Computation
% %%% semi-implicit Euler method %%%
% nSaveStep=floor(1/dt);
% nstep = 1;
% t = 0;
% figure(3)
% while (abs(u(n1))<thresh)
%     nstep = nstep+1;
%     u = U_U\(L_U\(u + dt*u.*g(u,x-c*t)));
%     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([-L L min(u)-0.1 max(abs(u))+0.1]) 
%         title(strcat('Time=' ,num2str(nstep*dt)),'fontsize',18)
%         drawnow;
%     end
%    t = t+dt;
% end
% 
% speed = [N*(Int1-Int2)*dx/(t_m-t_p) 2*sqrt(d*alpha)];
% 



