% This file solves the Fisher's equation
% on a square 2-dimensional rectangular grid with homogeneous 
%Dirichlet boundary conditions
 
clear all

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

d = 1; % diffusion coefficient = 1 wlog

Llist= [19,22]; %Length scale factor
N=30; % number of interior spatial grid points
% in each dimension 
h = 1/(N+1); %grid size
M = 18; % number of time steps to take

for j = 1:2 % simulate for each of the scale factors on the list
lam = Llist(j);
dt = h/5;
 scal = d*dt/h^2;
 X = h*(0:N+1)';
X0 = h*(1:N);
    cse = j;
    switch cse % specify initial data
        case 1   
        V0 =  ( tanh((X0' )/(1/15))).*( tanh(-(X0'- 1)/(1/15)));
        case 2
        V0 =  (exp(-(X0'-0.5).^2/(0.01)))/10;
 end
    
 figure((j-1)*(M+1)+1)  %plot initial data
  surf(X,X, [0;V0;0]*[0;V0;0]')
  axis([0 1 0 1 0 0.2])
  zlabel('v(\xi,\eta)','fontsize',20)
 V = reshape(V0 * V0',N^2,1)  ;
 Vmx(1) = max(V);
 
kstep = 50;
t=0;
 % uses Crank Nicolson for the diffusion part
% and explicit Euler for reaction part
Atm=(1+2*scal)*ones(N^2,1);  % the diagonal
A_off_d = ones(N-1,1);
for jp = 1:N-1
    A_off_d = [A_off_d;0;ones(N-1,1)];
end
Am = diag(Atm) - scal/2*(diag(ones(N*(N-1),1),N)+diag(ones(N*(N-1),1),-N)...
    +diag(A_off_d,-1)+diag(A_off_d,1));
 
% The algorithm
  
for n=2:M
for k = 1:kstep
 %F=V.*(1-V).*(V-0.1);  % bistable rhs
 %F =  V.*(1-V)-V.^2./(R*(ep^2+V.^2)); % spruce budworm rhs
 F =lam* V.*(1-V);  % Fisher equation rhs
  
 Vt = scal/2*([A_off_d;0].*[V(2:end);0] +[0;A_off_d].*[0;V(1:end-1)]...
     +[zeros(N,1);V(1:end-N)] + [V(N+1:end);zeros(N,1)]-4*V) +dt*F;
 
V = Am\(V+Vt);% update V
t=t + dt;
end

Vmx(n) = max(V);
T(n) = t;
figure((j-1)*(M+1)+n)
U = reshape(V,N,N);
Uh = [zeros(1,N+2);[zeros(N,1),U,zeros(N,1)];zeros(1,N+2)];
 surf(X,X,Uh)
 axis([0 1 0 1 0 0.2])
 
zlabel('v(\xi,\eta)','fontsize',20)
formatSpecF = '%6.2f\n';
title(strcat('Time = ',sprintf(formatSpecF,t)),'fontsize',18)
   

end
figure(j*(M+1))
plot(T,Vmx,'linewidth',2)
xlabel('time','fontsize',20)
ylabel('max(V)','fontsize',20)

end

