%  code to solve the two piece TBI problem
function tp_tbi
%First the phase portrait
global alp
%specify alpha
alp = 0.25;

N = 50;
du = 1/N;
u1 = [0:du:1];
w1 = sqrt(2*(F(1)-F(u1)));
w1a = sqrt(2*(F(1)-F(alp)));
Ar = 1-F(1)/F(alp);
umx = (2*alp)/3 + 2/3 - sqrt(4*alp^2 - 10*alp + 4)/3
du0=umx/N;
u0=[0:du0:umx];
w0=sqrt(-2*Ar*F(u0));

figure(1)
plot(u1,w1,u0,w0,alp,w1a,'k*')
xlabel('U','fontsize',20)
ylabel('W','fontsize',20)

%now make a plot of the critical curve
a=[00.01:0.01:0.5];
Ar = (1 + a).*(a - 1).^3./(a.^3.*(a - 2));

figure(2)
plot(a,Ar)
axis([0 0.5 0 100])
xlabel('\alpha','fontsize',20)
ylabel('A_r','fontsize',20)
text(0.25,60,'propagation blocked','fontsize',20)

function FI = F(u)
global alp

FI=-u.^4/4 + ((1 + alp)*u.^3)/3 - alp*u.^2/2;

