% This file solves equation u_t = u_xx + f(u)
 % using the Crank Nicholson method
 % Feel free to specify your own f in the function feval

 function CN_diff_rxn
global alf

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
 
 L = 30; % length of the domain
d = 1; % diffusion coefficient
alf = 0.25;
 Xc = L;  % center of initial profile
 %Xc = 0;
 %Xc = L/2;

 %to include a region of block
%uncomment the desired length of block region
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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Warning:  to get this simulation to work for your equation,
%the parameters N and dt may need to be adjusted.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N=300;   % number of spatial grid points
M = 61;  % number of time steps
h = L/N;
dt = h/2;
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 = 100;  % number of intermediate steps
t=0;
tn = zeros(M,1);
 Vs = zeros(N,M);
 
 a = 0.42;  % select a value of a, amplitude of the initial profile
 % a = 0.41;  % select a value of a, amplitude of the initial profile

a = 0.5;
lam = 4.2; % select a value of lam
  V = a*sech((X-Xc)/lam).^2; %initial profile

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

  
  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

% A 3d view of the wave
 figure(2)
 mesh(tn,X,Vs)
 ylabel('\xi', 'fontsize',20)
  xlabel('\tau', 'fontsize',20)
   zlabel('u(\xi,\tau)', 'fontsize',20)
   
   
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 the speed, this curve needs to be truncated
% some
spx = sfc(1, fix(sfcl/10):end-2);
spy = sfc(2, fix(sfcl/10):end-2);
 
  % 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)
  
  
     function F=feval(V)
     global alf
        
       F=V.*(1-V).*(V-alf);% bistable equation; 
       
 %F = V.*(1-V); %Fisher equation