% This file solves the bistable equation and determines threshold
% saves data to a file
function threshold
global alp

% because this code takes a while to execute, it is best to run it for
% different values of alpha,save the data, 
% then make plots using threshold_plotter.m at the end.

%
alplist = [0.1,0.25,0.4];
for ja = 1:3
alp =alplist(ja) 

for k = 1:50
    sc = k/5
    sck(k) = sc;
 
mpb = 0;
mpa = 1;

zoa = thr(sc,mpa);
if (zoa == 1)
    for j = 1:12
        mpc = (mpa+mpb)/2;
        zoc = thr(sc,mpc);
        mpa = (1-zoc)*mpa+zoc*mpc;
        mpb = zoc*mpb+(1-zoc)*mpc;
    end
   amp(k) =  mpc;
   mpc
else
    amp(k) = 1;
end

end
figure(1)
plot(sck,amp,'linewidth',2)
hold on
end

xlabel('length constant, \lambda','fontsize',20)
ylabel('Amplitude','fontsize',20)
text(5,0.8,'superthreshold','fontsize',20)
text(0.3,0.18,'subthreshold','fontsize',20)
 ja
switch ja
    case 1
save('thresh_0p1','sck', 'amp')
case 2
save('thresh_0p25','sck', 'amp')

case 3
save('thresh_0p4','sck', 'amp')


end
 


function zero_one=thr(sc,amp)
global alp
L = 100; % length of the domain
d = 1; % diffusion coefficient
 
N=200;   % number of spatial grid points
M = 39; 
h = L/N;
dt = h/2;
 
scal = d*dt/h^2;

X = h*(1:N)';
kstep = 10;
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;
%  
%  figure(1)
%  plot(X,V,'--')
%  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

end
 
zero_one = (V(1)>alp);
 