% This file solves equation u_t = u_xx + f(u)
 % using Crank Nicholson method
 % with f(u) from the spruce budworm model

 function CN_diff_rxn
global sigma kappa

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
 
% parameters
%
kappa = 1.11/355;

for js = 1:20 % 1:20
sigma = 5.4+(js-1)*8;
 
sig(js)= sigma;
 L = 35; % length of the domain
d = 1; % diffusion coefficient, nondimensional

Xc = 0;  % center of initial profile
% Xc = L/2;
% Xc = L;

%to include a region of block

%Lb = 2.95;  %size (length) of block region
%Lb = 3.0; %this length gives failure
Lb = 0.;  % there is no block region

Xb = 15; % location of block region
kd = 0.1;  % decay constant in block region
N=1000;   % number of spatial grid points
M = 30;  % number of time steps
h = L/N;
dt = h/10;
rat = ones(M,3);
 scal = d*dt/h^2;

X = h*(1:N)';
Xmsk = (X>=Xb)&(X<=(Xb+Lb));
Xblk = find((X>=Xb)&(X<=(Xb+Lb)));
t0 = L/(16*d);
 sc = 0.2*sqrt(4*pi*d*t0);
kstep = 10;  % number of intermediate steps
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);
a = 0.42;  % select a value of a, amplitude of the initial profile
a = 0.75;
lam = 2; % select a value of lam for the initial profile
  V = a*sech((X-L)/lam).^2+ 1.e-4;; %initial profile
   
  Vs(:,1) = V;
 figure(1)
 plot(X-Xb,V,'--','linewidth',2)
 hold on 
   
% uses Crank Nicolson to solve the diffusion equation
 % The algorithm
 
for n=2:M
for k = 1:kstep
   % your favorite function F goes here
   
  F = feval(V);
 Fb = Xmsk.*(-kd*V) + (1-Xmsk).*F;  % this accounts for the block region
V_t =  scal/2 * ([V(2:N );V(N )] - 2*V + [V(1 );V(1:N-1 )]) +dt*Fb ;
V  = Am\(V+V_t);
t=t + dt;

end

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

end
   
figure(1)
plot(X(Xblk)-Xb,0.01*ones(length(Xblk),1),'k','linewidth',5)
xlabel('\xi', 'fontsize',20)
  ylabel('u(\xi,\tau)', 'fontsize',20)
%  axis([-Xb L-Xb 0 1])
hold off

TL = t;
 figure(3)
sfc=contour(tn,X,Vs,[0.5 0.5],'linewidth',2);
sfcl=length(sfc(1,:));
%to get a good estimate of speed, need to truncate this curve some
spx = sfc(1, fix(sfcl/3):end-fix(sfcl/30));
spy = sfc(2, fix(sfcl/3):end-fix(sfcl/30));
 
 

  % now use sfc to estimate the speed
  % this must be done interactively
  
  p=polyfit(spx,spy,1) %p(1) is the estimated speed
  spest = p(2)+p(1)*spx;
  
  figure(3)
plot(spx,spy,spx,spest,'--')
 ylabel('\xi', 'fontsize',20)
  xlabel('\tau', 'fontsize',20)
  
  speed(js) = p(1);
end

figure(4)
plot(sig,-speed,'linewidth',2)
xlabel('\sigma','fontsize',24)
ylabel('speed','fontsize',24)
set(gca, 'XTick', 0:25: 150);
  
     function F=feval(V)
     global sigma kappa
        
       %F=V.*(1-V).*(V-alf);% bistable equation; 
     F = sigma*V.*(1-V)-V.^2./(kappa^2+V.^2);
     
 %F = V.*(1-V); %Fisher equation