%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); 
% 
% 
cmykflg=0;
%if cmykflg=1, this code will make cmyk figure files
%cmykflg = 1;

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

%guess a starting value for u0
 
a0list = [0.23, 0.2442, 0.25];
for j = 1:3
    a0 = a0list(j);
v = a0;
fa = v*(1-v)*(v-alf);

u0 = a0 -fa*t0^2/4; % use power series to find initial data
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)); 

% the area under the curve
N2 = 2*pi*S(end,3)
 
figure(1)
 
plot(S(:,1),S(:,2),'linewidth',2)
 
 xlabel('u','fontsize',20)
 ylabel('w','fontsize',20)
 axis([0 0.25 0 0.2])
 annotation('textarrow',[0.48,0.38 ],[0.85,0.85 ],'linewidth',2)
  hold on
  
 figure(2)
 plot(T,S(:,1),'linewidth',2)
 xlabel('\xi','fontsize',20)
 ylabel('u(\xi)','fontsize',20)
 axis([0 25 0 0.25])
 
hold on

end
figure(1)
plot([alf,alf],[0,0.2],'--')
legend('a = 0.23', 'a = 0.2442', 'a = 0.25','u = \alpha')

hold off
if(cmykflg==1) 
print('../../figs_c/chapt_8/two_d_thres_pp','-deps','-cmyk')
end


figure(2)
legend('a = 0.23', 'a = 0.2442', 'a = 0.25')
hold off
if(cmykflg==1) 
print('../../figs_c/chapt_8/two_d_thres_prf','-deps','-cmyk')
end


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


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