% calculate the threshold stimulus numbers
function bistable
global alf c
 
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
% 
alf = 0.1;
for j = 1:200
    alf = j/401;
    al(j) = alf;
 c = 0;
lam = (-c+sqrt(c^2+4*alf))/2;
u0 = 0.001;
w0 = lam*u0;

s = [u0;w0;0];
tstep = 0.1;
t_end = 250;

tspan = [0:tstep:t_end];
[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 

int(j) = 2*S(end,3); %doubled because we integrated only half the interval

 
end

figure(1)
plot(al,int,'linewidth',2)
xlabel('\alpha','fontsize',20)
ylabel('N_c','fontsize',20)
function s_prime=deRHS(t,s)
global alf c
v = s(1);
w = s(2);
 
Fv = w;
Fw = c*w-v*(1-v)*(v-alf);
Fu = v;  %quadrature

s_prime = [Fv Fw Fu]';


% event function
function [value,isterminal,direction]=events(t,s)
v = s(1);
w = s(2);

value = w*(1-v); 
isterminal = 1;
direction = 0;
