% This code solves the diffusion equation with quadratic reaction term
% (i.e., Fisher)

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

cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;


 L = 100; % length of the domain
  
d = 1; % diffusion coefficient  
kf = 1;  % time constant for the nonlinear term
 
N=150; % number of spatial grid points
M = 50; % total number of time steps
h = L/N;
 
 
kstep = 50; %number of intermediate time steps
dt = h/50; % the actual time step for CN

%dt = 0.005;
 scal = d*dt/h^2; % scale factor

X = h*(1:N)';  %the spatial grid
 t=0;
tn = zeros(M,1);
Vs = zeros(N,M);
% set up matrix

Atm=(1+scal)*ones(N,1);
Atm(1,1) = 1+scal/2;
Atm(N,1) = 1+scal/2;

Am=diag(Atm)-diag(scal/2*ones(N-1,1),1)-diag(scal/2*ones(N-1,1),-1);
 
 x0 = L/2;%
 %x0 = 0; %The center of the initial profile
  V = 0.002*exp(-(X-x0).^2/(L/4));  % initial data a small gaussian distribution
%V=exp(-0.05*X.^2);  % alternate initial data ; with these initial data
%increase number of time steps
 
  Vs(:,1) = V;
 figure(1)
 plot(X,V,'--','linewidth',2)
 hold on 
   
% uses Crank Nicolson to solve the diffusion equation
 % The algorithm:
for n=2:M
for k = 1:kstep
    
 % F = K*V.*(1-V).*(V-0.1);% right hand side for the bistable equation; 
 F = kf*V.*(1-V); %Fisher equation
 
V_t =  scal/2 * ([V(2:N );V(N )] - 2*V + [V(1 );V(1:N-1 )]) +dt*F ;
V  = Am\(V+V_t);
t=t + dt;

end

 figure(1)
 plot(X,V,'linewidth',2)
 Vs(:,n) = V;  % save V
  
 hold on
 pause(0.05)
tn(n) = t;  % save t

end
   
figure(1)
xlabel('x', 'fontsize',20)
  ylabel('u(x,t)', 'fontsize',20)
hold off
if(cmykflg==1) 
print('../../figs_c/chapt_6/fisher_t_dep','-deps','-cmyk')
end


 figure(2)
 mesh(tn,X,Vs)
 ylabel('x', 'fontsize',20)
  xlabel('t', 'fontsize',20)
   zlabel('u(x,t)', 'fontsize',20)
TL = t;
 

figure(3)
contour(tn,X,Vs,[0.5 0.5],'linewidth',2)  % make a contour plot of the V = 0.5 level surface
 ylabel('x', 'fontsize',20)
  xlabel('t', 'fontsize',20)
set(gca, 'XTick', 0 : 5 : 40); 

if(cmykflg==1) 
print('../../figs_c/chapt_6/fisher_xt_locus','-deps','-cmyk')
end
% now a little post processing of Vs to find the speed of travelling
  % waves
  if(x0==0)
   xf = zeros(M,1);
  for j = 1:M
      if (max(Vs(:,j)>=0.5))
      xj = max(find(Vs(:,j)>=0.5));
      % estimate the x value where Vs = 0.5
      if(xj<N)
   
      Xf(j) = X(xj) +(0.5-Vs(xj,j))*h/(Vs(xj+1,j)-Vs(xj,j));  %This estimates where V=0.5
      end
      end
      
  end
  
  jmin = min(find(Xf>0));
  tmin= tn(jmin);
  xmin = Xf(jmin);
  
  jmx = max(find(Xf<L));
  tmx= tn(jmx);
  xmx = Xf(jmx);
  
  rngcut = 15;
  rnge = [jmin+rngcut:jmx-rngcut];  % shorten the range a bit to get a better straight line
  
  p = polyfit(tn(rnge),Xf(rnge)',1)
  p(1) % is an estimate of the speed
  figure(4)
  plot(tn(rnge) ,Xf(rnge) ,tn(rnge),p(1)*tn(rnge)+p(2),'--','linewidth',2)
  end
  