%shooting to find the 3 d threshold profile
function three_d_threshold
global alf 

set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
% 
% 
alf = 0.1;
t0 = 1.e-6;  % a small starting value 

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

u0 = a0 -fa*t0^2/6; %Estimate initial data using  power series
w0 = fa*t0^2/3;

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

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

%this is the integral N_3
4*pi*S(end,3)

figure(1)
 
plot(S(:,1),S(:,2),'linewidth',2)
 
 xlabel('u','fontsize',20)
 ylabel('w','fontsize',20)
 axis([0 0.45 0 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.45])
 
hold on

end
figure(1)
plot([alf,alf],[0,2],'--')
legend('a = 0.37', 'a = 0.404', 'a = 0.41','u = \alpha')

hold off

figure(2)
legend('a = 0.37', 'a = 0.404', 'a = 0.41')
hold off



 
function s_prime=deRHS(t,s)
global alf 
u = s(1);
w = s(2);
 
Fu = -w/t^2;
Fw =  t^2*u*(1-u)*(u-alf);
Fn = t^2*u;  %the integral of the profile
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;
