% This file solves the bistable equation determines threshold


global alp

alp = 0.25;
sclist=[4.2,4.3];
amlist = [0.39,0.4];
for kk = 1:2
sc = sclist(kk);
amp = amlist(kk);
  
L = 50; % length of the domain
d = 1; % diffusion coefficient
 
N=101;   % number of spatial grid points
M = 19; 
h = L/(N-1);
dt = h/2;
 
scal = d*dt/h^2;

X = h*(0:N-1)';
kstep = 20;
t=0;
 
% 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);

% initial data
  V = amp./cosh(X/sc).^2;
  Vp = 2*amp*sinh(X/sc)./cosh(X/sc).^2/sc;
  figure(2)
  plot(V,Vp)
  
%  
 figure(kk)
 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
    
    
 % bistable equation; F = K*V.*(1-V).*(V-0.1);
 F = V.*(1-V).*(V-alp); %bistable 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(kk)
 plot(X,V)
 
 hold on
 pause(0.05)
%  

end
 
figure(kk)
axis([0 20 0 1])
xlabel('\xi', 'fontsize',20)
  ylabel('U(\xi,\tau)', 'fontsize',20)
hold off
 
end
