% barrier block for bistable equation
%This code calculates the function Y(U(0))
%this code has been modified for Wolbachia
function block
global alp  Ku kappa un F_int U2 F0
%specify alpha  
 
set(0,                           ...
   'defaultaxesfontsize', 20,   ...
   'defaultaxeslinewidth', 1.2, ...
   'defaultlinelinewidth', 2.0, ...
   'defaultpatchlinewidth', 0.7); 
 
alplist=[0.3]; % list of alpha values to work through
klist=[2]; % list of k values to work through
for kk = 1:length(klist)  %outer loop, go through the kappa values
kappa  = klist(kk);
 
ksq = sqrt(kappa);
for jk = 1:length(alplist)  %inner loop, go through all the alpha values
   
% pick a value of alpha
alp = alplist(jk);
aj(jk) = alp;

 %here are the guts of the algorithm:
%%Specify (or determine) the three zeros of f(u)
U0 = 0;
U1=alp;
U2 = 1;
N = 250; % discretize of the u interval
du = (U2-U0)/N;
un = [U0:du:U2];

% calculate the integral of f for all un
 F_int = cumsum(feval(un))*du;  % feval is where the bistable function is evaluated
 
 %find the zero of fint using bisection
 % guess a lower bound
 a = U0+1.e-6;
 % guess an upper bound
 b=U2;
 
umx = bisect(a,b,@fint)  %This is the zero of F_int
% which is  the maximum value for U(0):

%plot the two basic trajectories

%the upper integral curve
u1=un(2:end);
F1 = fint(U2);
w1 = sqrt(2* (F1-fint(u1)));
Fa = fint(U1);
 w1s=sqrt(2*(F1-Fa));  % the maximal value of the upper curve
% the lower integral  curve
du0 = (umx-1.e-6-U0)/N;
u0 = [U0:du0:umx];
Fw = fint(u0);
w0 = sqrt( -2*Fw); % these are starting values of w for given U(0)
w0s = sqrt(-2*fint(U1)); %the maximal value on the lower curve

Ku = w0.^2/2-kappa*u0.^2/2; % Ku is the integration constant
% for the connecting curve

 % use bisection to find the end value of u(Y))
 ua = u0;
 Fa = Fuy(ua);
 ub = ones(1,length(u0));
 Fb = Fuy(ub);
    uc = bisect(ua,ub,@Fuy);
 
% uc is the calculated end value of u(L)
% calculate L from the analytical solution for u
% since u satisfies u''-ku = 0, u(0) = u0, u'(0) = w0
A = (ksq*u0+w0)/(2*ksq);
B = (ksq*u0-w0)/(2*ksq);
yp = (uc+sqrt(uc.^2-4*B.*A))./(2*A);

Lp = log(yp)/ksq;

wc = ksq*A.*yp-ksq*B./yp;  %w(L)
  

[Lmin,jw] = min(Lp);
 
Lk(jk) = Lmin
  
     % you may want to pick a different value of alpha here
     
% for this particular value of L plot the profile
 % this will give an error if kappa is too small:
 
  y0 =   sqrt( abs(B(jw)./A(jw))); 
   
  x0=log(y0)/(ksq);
   
   dx = (Lmin-x0)/N;
   
x = [x0:dx:Lmin];
ux = A(jw)*exp(ksq*x)+B(jw)*exp(-ksq*x);
wx = A(jw)*ksq*exp(ksq*x)-ksq*B(jw)*exp(-ksq*x);
 
 t=[0:dx:Lmin];
 ut = A(jw)*exp(ksq*t)+B(jw)*exp(-ksq*t);  %This is the spatial trajectory 
 % in the block region, known analytically
 
figure(1)  % phase portrait of the blocking solution
  plot(u0,w0, u1,w1, ux,wx,'--',uc(jw),wc(jw),'k*',u0(jw),w0(jw),'k*',...
      'linewidth',2.5)
   text(0.45,0.057,'U(0)','fontsize',20)
   text(0.48,0.27,'U(Y)','fontsize',20)
   annotation('arrow',[.49 .5],[.43,.55])
   annotation('arrow',[.78 .82],[.4,.3])
   annotation('arrow',[.169 .24],[.16,.25])
  xlabel('U','fontsize',20)
   ylabel('W','fontsize',20)
   axis([0 1 0 0.3])
   
  
    figure(2)  %plot Y as a function of U(0)
    plot(u0,Lp,u0(jw),Lmin,'k*','linewidth',2)
    xlabel('U(0)','fontsize',20)
    ylabel('Y','fontsize',20)
    axis([0 umx 0 20])
 
 % for this value of alpha, create the spatial profile
% the upper trajectory
t_end = 7;
tstep = 0.1;

%specify the output points
tspan = [0:tstep:t_end];
 
   
s0=[uc(jw)];
F0=fint(U2);
[T1,S1] = ode23(@deRHS1,tspan, s0, odeset('maxstep',1)); 
% the lower trajectory
 t_end = -10;
tstep = -0.1;

%specify the output points
tspan = [0:tstep:t_end];

 s0 = [u0(jw)];
 F0 = fint(U0);
[T0,S0] = ode23(@deRHS1,tspan, s0, odeset('maxstep',1));  

figure(5) 

plot(T0,S0(:,1),T1+Lmin,S1(:,1),t,ut,'--',...
    0,u0(jw),'k*',Lmin,uc(jw),'k*','linewidth',2)
xlabel('\xi','fontsize',20)
ylabel('U','fontsize',20)
text(0.39,0.4,'(0,U(0))','fontsize',20)
text(0.45,0.5,'(Y,U(Y))','fontsize',20)
axis([-10 10 0 1])
  end 
end

 
 


function s_prime=deRHS(t,s)  % right hand side of the differntial equation
 
 u = s(1);
 w = s(2);
 Fu = w;
 Fw = -feval(u);
 s_prime = [Fu;Fw];
 
 function s_prime=deRHS1(t,s)  % right hand side of the differntial equation
 global F0
 u = s(1);

 Fu = sqrt(2* (F0-fint(u)));
 
 s_prime = [Fu];

 % here is where to specify the function f(u);
function f=feval(u);
global alp
f = u.*(1-u).*(u-alp);  % whatever function you want goes here

        
function fI = fint(u)
global un F_int

% here we assume that no analytical formula is available
%Calculate the integral of f using the vector F_int calculated at first
% u may be a vector
for j = 1:length(u)
ndx = min(find(u(j)<=un));
if(ndx==1)
    fI(j) = 0;
else
%linear interpolation
u1 = un(ndx-1);
f1 = F_int(ndx-1);
u2 = un(ndx);
f2 = F_int(ndx);
fI(j) = ((u(j)-u1)*f2-(u(j)-u2)*f1)/(u2-u1);
end
end

function fcross = Fuy(u)
global kappa Ku U2
fcross = kappa*u.^2/2+fint(u)-fint(U2)+Ku;

function root = bisect(a,b,feval)

ul = a;
fl = feval(ul);
uu = b;
fu = feval(uu);

% we make the assumption, without checking, that 
% fu*fl<0

% if not, the algoritm fails to find a root.

N = 20;  % number of iterates
% the main bisection algorithm
for j = 1:N
uc = (ul+uu)/2;
fc = feval(uc);
ftest = (fc.*fl>=0);
ul = ftest.*uc+(1-ftest).*ul;
fl = ftest.*fc + (1-ftest).*fl;

uu = (1-ftest).*uc+ ftest.*uu;
fu = (1-ftest).*fc + ftest.*fu;
end
root = uc;
