%shooting to find the 2 d threshold profile
function two_d_threshold
global alf 

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 

eps = 1.e-6;
% 
alf = 0.3;
t0 = eps;  % a small starting value 

%guess a starting value for u0
amin = alf + eps
amx = 1-eps
flag = 0;
 a0 = (amin+amx)/2;
for j = 1:20
   
     
ans = homclin(a0,flag)
if(ans == 1) amx = a0;
end
if (ans == 2) amin = a0;
end
a0 = (amin+amx)/2

%ans = 1 means a0 too high
%ans = 2 means a0 too low

end

%
flag = 1;
ans = homclin(a0,flag)

function shoot = homclin(a0,flag)
global alf

t0 = 1.e-6;  % a small starting value 

%guess a starting value for u0
v = a0;
fa = v*(1-v)*(v-alf);

u0 = a0 -fa*t0^2/4;
w0 = -fa*t0^2/2;

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

tspan = [ t0:tstep:t_end];
[T,S] = ode23s(@deRHS,tspan,s,odeset('Events',@events)); 
 if (flag ==1)
     figure(1)
     plot(T,S(:,1),'linewidth',2)
     xlabel('\xi','fontsize',20)
     ylabel('u(\xi)','fontsize',20)
     figure(2)
     plot(S(:,1),-S(:,2),'linewidth',2)
     xlabel('u','fontsize',20)
      ylabel('w','fontsize',20)
     N2 = 2*pi*S(end,3)
 end
 
 [A,j] = min([abs(S(end,1)),abs(S(end,2))]);
 shoot = j;

 
function s_prime=deRHS(t,s)
global alf 
u = s(1);
w = s(2);
 
Fu = w/t;
Fw =  -t*u*(1-u)*(u-alf);
Fint = t*u; 
s_prime = [Fu Fw Fint]';


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

% integrate until one of w or v changes sign
value = w*v; 
isterminal = 1;
direction = 0;
