% This file solves the Fisher's equation
% on a two-dimensional circular domain with homogeneous 
%Dirichlet boundary conditions with radial symmetry


clear
Llist = [2;5];
d = 1; % diffusion coefficient = 1 wlog
 
N=91; % number of radial grid points
M = 40; % number of time steps

for j = 1:2
    
 L = Llist(j);
 h = 1/(N+1);
dt = h/5;
 scal = d*dt/h^2;
 
 X = h*(0:N+1)';  %the radial variable
X0 = h*(1:N);  %the interior grid points
    cse = j;
    switch cse % specify initial data
        case 1
       V =  1 -tanh((X0'-0.8)/(1/20)) ;
        case 2
        V = 0.01* (1-tanh((X0.^2' )/(1/20))) ;
    end

 figure(j)
  plot(L*X,[V(1);V;0],'--','linewidth',3)
   
  hold on
kstep = 10;
t=0;
 
% now make the radially symmetric laplacian
% uses Crank Nicolson for the diffusion part
% and explicit Euler for reaction part
Atm=(1+scal)*ones(N,1);  % the diagonal entry
Atm(1,1) = 1+3*scal/4; % this gives zero derivative at the origin
  
%  this is for the radial laplacian operator
Am=diag(Atm)-diag(scal/2*([1:N-1]+1/2)./[1:N-1],1)...
    -diag(scal/2*([1:N-1]+1/2)./[2:N],-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 = L^2*V.*(1-V);  % Fisher's equation
 
Vt =  scal/2 * ((([1:N]+1/2)./[1:N])'.*[V(2:N );0]...
    - 2*V + (([1:N]-1/2)./[1:N])'.*[V(1);V(1:N-1)]) +dt*F  ;
 
V = Am\(V+Vt);
t=t + dt;
 
end
figure(j)
  plot(L*X,[V(1);V;0],'linewidth',2)
  %axis([0 L  0 1])
 
xlabel('r','fontsize',20)
ylabel('v(r,\tau)','fontsize',20)
   hold on

end

end

