function [B0]=intvalWrightCont(B0)

% The input is B0=[alpha0;r0;x0]

count=0;   % Total Number of Continuation Step(s) 

s=3;
alpha0=B0(1,1);
dal=5e-5;
dal_min=1e-15;
dal_max=2;
dal0=0;

hx=B0(3:length(B0));
m=length(hx)/2;
B0_intermediate=B0;

change_dimension=0;

%%%%%%%%%%%%%%%%%%%%%%
% Start continuation %
%%%%%%%%%%%%%%%%%%%%%%

while count<300 %(alpha0<alpha_max),  %%% until we reach the endpoint of the continuation %%%
    
    if change_dimension>0
        B0_intermediate((2*m+3:2*m+4),:)=0;
        display('We increased the dimension by 2')
        m=m+1;
        hx=B0_intermediate((3:2*m+2),count+1);
        change_dimension=0;
    end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                          %%%
%%% (Quasi-Newton iteration) %%%
%%%                          %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
tol=1.e-14; 
x1=hx;  
F=F_Wright(x1,alpha0);
J=inv(F_Wright_x(x1,alpha0));
itercount=0;
while (itercount<=7) && (norm(F) > tol),
    x1=x1-J*F;                   %%% Quasi-Newton Corrector %%%
    F=F_Wright(x1,alpha0);
    J=inv(F_Wright_x(x1,alpha0));
    itercount=itercount+1;                            
end
if norm(F)>tol                   %%% The corrector did not converged %%%
  display(['Failure! Problem with Quasi-Newton predictor-corrector at alpha = ',int2str(alpha0)])
  return
end  

bx=x1;
dx=-J*F_Wright_alpha(bx);    %%%  Tangent Vector  %%%   

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                       %%%
%%% CONSTRUCTION OF THE RADII POLYNOMIALS %%%
%%%                                       %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       
i0=intval('0');
i1=intval('1');
i2=intval('2');
i3=intval('3');
i4=intval('4');
i6=intval('6');
i8=intval('8');
ipi=midrad(3.14159265358979323,1e-16);
ialpha0=intval(alpha0);

bL=intval(bx(1,1));
ba=intval([bx(2,1);bx((2*((1:m-1))'+1),1)]);
bb=intval([0;bx((2*((2:m))'),1)]);
dL=intval(dx(1,1));
da=intval([dx(2,1);dx((2*((1:m-1))'+1),1)]);
db=intval([0;dx((2*((2:m))'),1)]);

Df_at_xF=F_Wright_x(intval([bx;zeros(2*m-2,1)]),ialpha0);
J_F=inv(Df_at_xF);

%%%%%%%%%%%%%%%%%%%%%
%%%               %%%
%%% D, C & hat{C} %%%
%%%               %%%
%%%%%%%%%%%%%%%%%%%%%

d0=F_Wright(intval([bx;zeros(2*m-2,1)]),ialpha0);
d1=Df_at_xF*intval([dx;zeros(2*m-2,1)])+F_Wright_alpha(intval([bx;zeros(2*m-2,1)]));

% Preallocate memory %
D2=intval(zeros(2*(2*m-1),1));
D3=intval(zeros(2*(2*m-1),1));
D4=intval(zeros(2*(2*m-1),1));
D5=intval(zeros(2*(2*m-1),1));
C10=intval(zeros(2*(2*m-1),1));
C11=intval(zeros(2*(2*m-1),1));
C12=intval(zeros(2*(2*m-1),1));
C13=intval(zeros(2*(2*m-1),1));
C20=intval(zeros(2*(2*m-1),1));
C21=intval(zeros(2*(2*m-1),1));
C22=intval(zeros(2*(2*m-1),1));
C30=intval(zeros(2*(2*m-1),1));
C31=intval(zeros(2*(2*m-1),1));

v_Ftemp=intval(zeros(4*(m-1),1));
v_Ftemp((2*(1:2*m-2)),:)=intval(1:2*m-2)'.^(-s);
v_Ftemp((2*(1:2*m-2))-1,:)=intval(1:2*m-2)'.^(-s);
v_F=[i1;i1;v_Ftemp];

C1_minus1=i2/((s-i1)*(i2*m-i2)^(s-i1));

for k=0:2*m-2
    
    sC10=0; 
    sD2=0; 
    sD3=0; 
    sD4=0; 
    sD5=0; 
    sC11_1=0; 
    sC11_2=0;  
    sC12_1=0; 
    sC12_2=0; 
    sC13=0;
    sC20_1=0; 
    sC20_2=0; 
    sC21=0;
    sC22=0;
   
    for k1=2*m-1:k+m-1
        sC10=sC10+i4*ialpha0*(abs(a(ba,k-k1))+abs(b(bb,k-k1)))/omega(k1,s);
    end
    
    for k1=-m+1:m-1
        
        if abs(k-k1)<m
           
        tsD2=abs(-(-(i1/i2)*k1^2*dL^2*(a(ba,k1)*a(ba,k-k1)-b(bb,k1)*b(bb,k-k1))+a(da,k1)*a(da,k-k1)-b(db,k1)*b(db,k-k1)+k1*dL*(a(ba,k1)*b(db,k-k1)+a(da,k1)*b(bb,k-k1)+b(bb,k1)*a(da,k-k1)+b(db,k1)*a(ba,k-k1)))*ialpha0-k1*dL*(a(ba,k1)*b(bb,k-k1)+b(bb,k1)*a(ba,k-k1))-a(ba,k1)*a(da,k-k1)-a(da,k1)*a(ba,k-k1)+b(bb,k1)*b(db,k-k1)+b(db,k1)*b(bb,k-k1))+abs(-(i1/i2)*ialpha0*(-k1^2*dL^2*(a(ba,k1)*b(bb,k-k1)+b(bb,k1)*a(ba,k-k1))-i2*k1*dL*(a(ba,k1)*a(da,k-k1)+a(da,k1)*a(ba,k-k1)-b(bb,k1)*b(db,k-k1)-b(db,k1)*b(bb,k-k1))+i2*a(da,k1)*b(db,k-k1)+i2*b(db,k1)*a(da,k-k1))+k1*dL*(a(ba,k1)*a(ba,k-k1)-b(bb,k1)*b(bb,k-k1))-a(ba,k1)*b(db,k-k1)-a(da,k1)*b(bb,k-k1)-b(bb,k1)*a(da,k-k1)-b(db,k1)*a(ba,k-k1));
        tsD3=abs((i1/i2)*k1^2*dL^2*(a(ba,k1)*a(ba,k-k1)-b(bb,k1)*b(bb,k-k1))-i3*a(da,k1)*a(da,k-k1)+i3*b(db,k1)*b(db,k-k1)-i2*k1*dL*(a(ba,k1)*b(db,k-k1)+a(da,k1)*b(bb,k-k1)+b(bb,k1)*a(da,k-k1)+b(db,k1)*a(ba,k-k1))-(i1/i2)*ialpha0*(i2*k1*dL*(i2*a(da,k1)*b(db,k-k1)+i2*b(db,k1)*a(da,k-k1))-k1^2*dL^2*(a(ba,k1)*a(da,k-k1)+a(da,k1)*a(ba,k-k1)-b(bb,k1)*b(db,k-k1)-b(db,k1)*b(bb,k-k1))))+abs((k1*dL*(i2*a(da,k1)*a(da,k-k1)-i2*b(db,k1)*b(db,k-k1))+(i1/i2)*k1^2*dL^2*(a(ba,k1)*b(db,k-k1)+a(da,k1)*b(bb,k-k1)+b(bb,k1)*a(da,k-k1)+b(db,k1)*a(ba,k-k1)))*ialpha0+i2*k1*dL*(a(ba,k1)*a(da,k-k1)+a(da,k1)*a(ba,k-k1)-b(bb,k1)*b(db,k-k1)-b(db,k1)*b(bb,k-k1))+(i1/i2)*k1^2*dL^2*(a(ba,k1)*b(bb,k-k1)+b(bb,k1)*a(ba,k-k1))-i3*a(da,k1)*b(db,k-k1)-i3*b(db,k1)*a(da,k-k1));
        tsD4=abs((i1/i2)*k1^2*dL^2*(a(ba,k1)*a(da,k-k1)+a(da,k1)*a(ba,k-k1)-b(bb,k1)*b(db,k-k1)-b(db,k1)*b(bb,k-k1))-(i1/i2)*ialpha0*(-a(da,k1)*a(da,k-k1)+b(db,k1)*b(db,k-k1))*k1^2*dL^2-k1*dL*(i2*a(da,k1)*b(db,k-k1)+i2*b(db,k1)*a(da,k-k1))+(-a(da,k1)*b(db,k-k1)-b(db,k1)*a(da,k-k1))*k1*dL)+abs((i1/i2)*k1^2*dL^2*(a(ba,k1)*b(db,k-k1)+a(da,k1)*b(bb,k-k1)+b(bb,k1)*a(da,k-k1)+b(db,k1)*a(ba,k-k1))-(i1/i2)*ialpha0*(-a(da,k1)*b(db,k-k1)-b(db,k1)*a(da,k-k1))*k1^2*dL^2+k1*dL*(i2*a(da,k1)*a(da,k-k1)-i2*b(db,k1)*b(db,k-k1))-(-a(da,k1)*a(da,k-k1)+b(db,k1)*b(db,k-k1))*k1*dL);
        tsD5=(i1/i2)*abs((-a(da,k1)*a(da,k-k1)+b(db,k1)*b(db,k-k1))*k1^2*dL^2)+(i1/i2)*abs((-a(da,k1)*b(db,k-k1)-b(db,k1)*a(da,k-k1))*k1^2*dL^2);

        tsC11_1=abs(k1*(a(ba,k1)*b(bb,k-k1)+b(bb,k1)*a(ba,k-k1)+ialpha0*(a(ba,k1)*b(db,k-k1)+a(da,k1)*b(bb,k-k1)+b(bb,k1)*a(da,k-k1)+b(db,k1)*a(ba,k-k1))))+abs(k1*(a(ba,k1)*a(ba,k-k1)-b(bb,k1)*b(bb,k-k1)+ialpha0*(a(ba,k1)*a(da,k-k1)+a(da,k1)*a(ba,k-k1)-b(bb,k1)*b(db,k-k1)-b(db,k1)*b(bb,k-k1))))+ialpha0*k1^2*abs(dL)*(abs(a(ba,k1)*a(ba,k-k1)-b(bb,k1)*b(bb,k-k1))+abs(a(ba,k1)*b(bb,k-k1)+b(bb,k1)*a(ba,k-k1)))+i2*ialpha0*abs(dL)*(abs(k1)+abs(k-k1))*(abs(a(ba,k1))+abs(b(bb,k1)))/omega(k-k1,s);
        tsC12_1=abs(k1*(a(ba,k1)*b(db,k-k1)+a(da,k1)*b(bb,k-k1)+b(bb,k1)*a(da,k-k1)+b(db,k1)*a(ba,k-k1)+ialpha0*(a(da,k1)*b(db,k-k1)+b(db,k1)*a(da,k-k1))))+abs(k1*(a(ba,k1)*a(da,k-k1)+a(da,k1)*a(ba,k-k1)-b(bb,k1)*b(db,k-k1)-b(db,k1)*b(bb,k-k1)+ialpha0*(a(da,k1)*a(da,k-k1)-b(db,k1)*b(db,k-k1))));
        tsC13=abs(k1)*(abs(a(da,k1)*a(da,k-k1)-b(db,k1)*b(db,k-k1))+abs(a(da,k1)*b(db,k-k1)+b(db,k1)*a(da,k-k1)));
        tsC20_1=ialpha0*k1^2*(abs(a(ba,k1)*a(ba,k-k1)-b(bb,k1)*b(bb,k-k1))+abs(a(ba,k1)*b(bb,k-k1)+b(bb,k1)*a(ba,k-k1)))+i2*ialpha0*(abs(k1)+abs(k-k1))*(abs(a(ba,k1))+abs(b(bb,k1)))/omega(k-k1,s);
        
        sD2=sD2+tsD2;
        sD3=sD3+tsD3;
        sD4=sD4+tsD4;
        sD5=sD5+tsD5;
        
        sC11_1=sC11_1+tsC11_1;
        sC12_1=sC12_1+tsC12_1;
        sC13=sC13+tsC13;
        sC20_1=sC20_1+tsC20_1;
            
        end
        
        tsC11_2=(i4*(abs(a(ba,k1)+ialpha0*a(da,k1))+abs(b(bb,k1)+ialpha0*b(db,k1))))/omega(k-k1,s);
        tsC12_2=(i4*abs(a(da,k1))+i4*abs(b(db,k1)))/omega(k-k1,s);
        tsC20_2=i2*ialpha0*(abs(k1)+abs(k-k1))*(abs(a(ba,k1))+abs(b(bb,k1)))/omega(k-k1,s);
        tsC21=(i2*(abs(k1)+abs(k-k1)))*(abs(a(ba,k1)+ialpha0*a(da,k1))+abs(b(bb,k1)+ialpha0*b(db,k1)))/omega(k-k1,s);
        tsC22=(i2*(abs(k1)+abs(k-k1)))*(abs(a(da,k1))+abs(b(db,k1)))/omega(k-k1,s);

        sC11_2=sC11_2+tsC11_2;
        sC12_2=sC12_2+tsC12_2;
        sC20_2=sC20_2+tsC20_2;
        sC21=sC21+tsC21;
        sC22=sC22+tsC22;
        
    end
    
    tD2=k*abs(dL)*(abs(a(da,k))+abs(b(db,k)))+abs(-a(da,k)+(i1/i2)*ialpha0*k^2*dL^2*a(ba,k)-k*dL*b(bb,k)-ialpha0*k*dL*b(db,k))+abs(k*dL*a(ba,k)+(i1/i2)*ialpha0*k^2*dL^2*b(bb,k)+ialpha0*k*dL*a(da,k)-b(db,k));
    tD3=abs(-i2*k*dL*b(db,k)+(i1/i2)*ialpha0*k^2*dL^2*a(da,k)+(i1/i2)*k^2*dL^2*a(ba,k))+abs(2*k*dL*a(da,k)+(i1/i2)*ialpha0*k^2*dL^2*b(db,k)+(i1/i2)*k^2*dL^2*b(bb,k));
    tD4=(i1/i2)*abs(k^2*dL^2*a(da,k))+(i1/i2)*abs(k^2*dL^2*b(db,k));
    
    D2(2*k+1,1)=tD2+sD2;
    D2(2*k+2,1)=D2(2*k+1,1);
    D3(2*k+1,1)=tD3+sD3;
    D3(2*k+2,1)=D3(2*k+1,1);
    D4(2*k+1,1)=tD4+sD4;
    D4(2*k+2,1)=D4(2*k+1,1);
    D5(2*k+1,1)=sD5;
    D5(2*k+2,1)=D5(2*k+1,1);  
    
    phi=sum((intval(1:k-1)'.*intval(k-(1:k-1)')).^(-s));
    SUMwk1wk2=phi+i1/omega(k,s)*(i4+i2/(s-i1));
    SUMk1wk1wk2=k*phi/i2+(k+i1)^(-1)*(i1+(i2*s-i3)^(-1))+k/omega(k,s)+(k+i1)^(-s+i1)*(i1+(s-i1)^(-1));
    
    tC11=k*abs(ialpha0*b(db,k)+b(bb,k))+k*abs(ialpha0*a(da,k)+a(ba,k))+i2/omega(k,s)+ialpha0*k^2*abs(dL)*(abs(a(ba,k))+abs(b(bb,k)))+i2*ialpha0*k*abs(dL)/omega(k,s)+k*abs(b(db,k))+k*abs(dL)/omega(k,s);
    tC12=k*(abs(b(db,k))+abs(a(da,k)));
    tC20=i4*ialpha0*k/omega(k,s)+i8*ialpha0*SUMwk1wk2+ialpha0*k^2*(abs(a(ba,k))+abs(b(bb,k)))+i2*k/omega(k,s);
    tC21=i2*k/omega(k,s)+i8*SUMwk1wk2;
    tC30=i4*ialpha0*SUMk1wk1wk2;
    tC31=i4*SUMk1wk1wk2;
    
    C10(2*k+1,1)=sC10;
    C10(2*k+2,1)=C10(2*k+1,1);
    C11(2*k+1,1)=tC11+sC11_1+sC11_2;
    C11(2*k+2,1)=C11(2*k+1,1);
    C12(2*k+1,1)=tC12+sC12_1+sC12_2;
    C12(2*k+2,1)=C12(2*k+1,1);
    C13(2*k+1,1)=sC13;
    C13(2*k+2,1)=C13(2*k+1,1);
    C20(2*k+1,1)=tC20+sC20_1+sC20_2;
    C20(2*k+2,1)=C20(2*k+1,1);
    C21(2*k+1,1)=tC21+sC21;
    C21(2*k+2,1)=C21(2*k+1,1);
    C22(2*k+1,1)=sC22;
    C22(2*k+2,1)=C22(2*k+1,1);
    C30(2*k+1,1)=tC30;
    C30(2*k+2,1)=C30(2*k+1,1);
    C31(2*k+1,1)=tC31;
    C31(2*k+2,1)=C31(2*k+1,1);
    
end

C10(1,1)=C1_minus1;
C11(1,1)=0;
C12(1,1)=0;
C13(1,1)=0;
C20(1,1)=0;
C21(1,1)=0;
C22(1,1)=0;
C30(1,1)=0;
C31(1,1)=0;

D2(1,1)=0;
D3(1,1)=0;
D4(1,1)=0;
D5(1,1)=0;

Y0=abs(J_F*d0);
Y1=abs(J_F*d1);
Y2=abs(J_F)*D2;
Y3=abs(J_F)*D3;
Y4=abs(J_F)*D4;
Y5=abs(J_F)*D5;

temp_matrix=intval(eye(2*(2*m-1)))-J_F*Df_at_xF;
if sup(norm(temp_matrix,inf))>=1
    display('J_M is not invertible!')
    %B0=B0_intermediate(:,count);
    return
end
I_JDf=abs(temp_matrix)*v_F;

inv_condition=(ialpha0*abs(i1+a(ba,0))/bL+i1)/i2;
if m <= sup(inv_condition)
   display('A is not invertible!')
   return
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construction of the hatC^{(l1,l2)} %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gamma=i2*((i2*m-i1)/(i2*m-i2))^s+(i4*log(i2*m-i3)/(i2*m-i1)+(ipi^2-i6)/i3)*(i2/(i2*m-i1)+i1/i2)^(s-i2);

shatC10=0; 
shatC11=0; 
shatC12=0; 
shatC20=0; 
shatC21=0;  
shatC22=0; 

for k1=1:m-1
    
    tshatC10=(i4*ialpha0/(i2*m-i1))*(abs(a(ba,k1))+abs(b(bb,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^s));
    tshatC11=(i4/(i2*m-i1))*(abs(a(ba,k1)+ialpha0*a(da,k1))+abs(b(bb,k1)+ialpha0*b(db,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^s));
    tshatC12=(i4/(i2*m-i1))*(abs(a(da,k1))+abs(b(db,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^s));
    tshatC20=(i2*ialpha0*k1/(i2*m-i1))*(abs(a(ba,k1))+abs(b(bb,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^s))+(i2*abs(a(ba,k1))+i2*abs(b(bb,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^(s-i1)))*ialpha0;
    tshatC21=(i2*k1/(i2*m-i1))*(abs(a(ba,k1)+ialpha0*a(da,k1))+abs(b(bb,k1)+ialpha0*b(db,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^s))+(i2*abs(a(ba,k1)+ialpha0*a(da,k1))+i2*abs(b(bb,k1)+ialpha0*b(db,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^(s-i1)));
    tshatC22=(i2*k1/(i2*m-i1))*(abs(a(da,k1))+abs(b(db,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^s))+(i2*abs(a(ba,k1))+i2*abs(b(bb,k1)))*(i1+i1/(i1-(k1/(i2*m-i1))^(s-i1)));
    
    shatC10=shatC10+tshatC10;
    shatC11=shatC11+tshatC11;
    shatC12=shatC12+tshatC12;
    shatC20=shatC20+tshatC20;
    shatC21=shatC21+tshatC21;
    shatC22=shatC22+tshatC22;
    
end

thatC11=i2/(i2*m-i1)+i2*ialpha0*abs(dL)*abs(i1+a(ba,0))+abs(dL)+i4*abs(a(ba,0)+ialpha0*a(da,0))/(i2*m-i1);
thatC12=i4*abs(a(da,0))/(i2*m-i1);
thatC20=i2*ialpha0+i4*ialpha0*abs(i1+a(ba,0))+(i8*ialpha0/(i2*m-i1))*(i4+i2/(s-i1)+gamma)+i2*ialpha0*abs(a(ba,0));
thatC21=i2+(i8/(i2*m-i1))*(i4+i2/(s-i1)+gamma)+i2*abs(a(ba,0)+ialpha0*a(da,0));
thatC22=i2*abs(a(da,0));
thatC30=i4*ialpha0*(i3+i2/(s-i1)+(i1/i2)*gamma);
thatC31=i4*(i3+i2/(s-i1)+(i1/i2)*gamma);

hatC10=shatC10;
hatC11=thatC11+shatC11;
hatC12=thatC12+shatC12;
hatC13=i0;
hatC20=thatC20+shatC20;
hatC21=thatC21+shatC21;
hatC22=thatC22+shatC22;
hatC30=thatC30;
hatC31=thatC31;

rho=(i2*m-i1)/((i2*m-i1)*bL-ialpha0*abs(i1+ba(1,1)));

p=intval(zeros(2*(2*m-1)+1,4));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                              % 
% We first verify if we can solve the polynomials with dal=0,  %
% in which case we compute B_1(r_min) and B_1(r_max) and prove %
% that the branch continues.                                   %
%                                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Y=Y0;

Z=intval(zeros(2*(2*m-1),3)); % Preallocate Memory 

Z(:,1)=I_JDf+abs(J_F)*C10;
Z(:,2)=abs(J_F)*C20;
Z(:,3)=abs(J_F)*C30;

hatZ(1,1)=i1/(i2*m-i1)^s*(rho^2*ialpha0/(i2*m-i1)*(abs(ba(1,1))+abs(i1+ba(1,1)))+rho)*hatC10;
hatZ(1,2)=i1/(i2*m-i1)^s*(rho^2*ialpha0/(i2*m-i1)*(abs(ba(1,1))+abs(i1+ba(1,1)))+rho)*hatC20;
hatZ(1,3)=i1/(i2*m-i1)^s*(rho^2*ialpha0/(i2*m-i1)*(abs(ba(1,1))+abs(i1+ba(1,1)))+rho)*hatC30;

p(1:2*(2*m-1),1:4) = [Y Z];
p(2*(2*m-1)+1,1:4) = [i0 hatZ];

p(:,2)=p(:,2)-[v_F;i1/(i2*m-i1)^s];   

I=[0,Inf];
k=-1;
while k<=2*(2*m-1)-1 && I(2)>I(1), 
  I_k=rootsfinder_radii_polynomial(mid(p(k+2,:)));
  I(1)=max(I(1),I_k(1));
  I(2)=min(I(2),I_k(2));
  k=k+1;
end

if I(1)<I(2)       %%% Check that radii polynomials are really negative %%%
   r1_min=intval(min(3*I(1),(I(1)+I(2))/2));
   r1_max=intval(max((I(1)+I(2))/2,0.5*I(2)));
   pr_min=p(:,4);
   pr_max=p(:,4);
   for j=3:-1:1
     pr_min=p(:,j)+r1_min*pr_min;
     pr_max=p(:,j)+r1_max*pr_max;
   end
   if max(sup(pr_min))>=0 || max(sup(pr_max))>=0
      display(['Failure! Problem: max(pr)>=0 with dal=0 at alpha = ',num2str(alpha0)])
      B0=B0_intermediate(:,count);
      return
   end
else
    display(['Failure! Problem with dal=0 at alpha = ',num2str(alpha0)])
    B0=B0_intermediate(:,count);
    return
end

%%% Rigorous branch following verification: Do C0 and C1 connect ? %%%
if count>0
r_B0=intval(B0_intermediate(2,1));
continuation_check1(1,1)=r1_max-r_B0; % Comparison of the tails of B0 and B1(r1_max)
continuation_check1(2:2*m+1,1)=intval(bx)+r1_max*v_F(1:2*m)-(intval(hx)+r_B0*v_F(1:2*m));
continuation_check1(2*m+2:4*m+1,1)=intval(hx)-r_B0*v_F(1:2*m)-(intval(bx)-r1_max*v_F(1:2*m));
if min(inf(continuation_check1)) < 0
    continuation_check2(1,1)=r_B0-r1_min; % Comparison of the tails of B0 and B1(r1_min)
    continuation_check2(2:2*m+1,1)=intval(hx)+r_B0*v_F(1:2*m)-(intval(bx)+r1_min*v_F(1:2*m));
    continuation_check2(2*m+2:4*m+1,1)=intval(bx)-r1_min*v_F(1:2*m)-(intval(hx)-r_B0*v_F(1:2*m));
    if min(inf(continuation_check2)) < 0
        display(['Failure! Problem in proving that C0 and C1 connect at alpha = ',num2str(alpha0)])
        return
    end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                               %
% Now that we know that the polynomials are solved for dal=0,   %
% we want to maximize the dal>0 for which they are still solved %
%                                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NEXTSTEP=0; 
while NEXTSTEP==0,     %%% As long as NEXTSTEP=0, you keep going %%%

Y=Y0+intval(dal)*Y1+intval(dal)^2*Y2+intval(dal)^3*Y3+intval(dal)^4*Y4+intval(dal)^5*Y5;

Z(:,1)=I_JDf+abs(J_F)*(C10+intval(dal)*C11+intval(dal)^2*C12+intval(dal)^3*C13);
Z(:,2)=abs(J_F)*(C20+intval(dal)*C21+intval(dal)^2*C22);
Z(:,3)=abs(J_F)*(C30+intval(dal)*C31);

hatZ(1,1)=i1/(i2*m-i1)^s*(rho^2*alpha0/(i2*m-i1)*(abs(ba(1,1))+abs(i1+ba(1,1)))+rho)*(hatC10+intval(dal)*hatC11+intval(dal)^2*hatC12+intval(dal)^3*hatC13);
hatZ(1,2)=i1/(i2*m-i1)^s*(rho^2*alpha0/(i2*m-i1)*(abs(ba(1,1))+abs(i1+ba(1,1)))+rho)*(hatC20+intval(dal)*hatC21+intval(dal)^2*hatC22);
hatZ(1,3)=i1/(i2*m-i1)^s*(rho^2*alpha0/(i2*m-i1)*(abs(ba(1,1))+abs(i1+ba(1,1)))+rho)*(hatC30+intval(dal)*hatC31);

p(1:2*(2*m-1),1:4) = [Y Z];
p(2*(2*m-1)+1,1:4) = [i0 hatZ];

p(:,2)=p(:,2)-[v_F;i1/(i2*m-i1)^s];   

I=[0,Inf];
k=-1;
while I(2)>I(1) && k<=2*(2*m-1)-1,
  I_k=rootsfinder_radii_polynomial(mid(p(k+2,:)));
  I(1)=max(I(1),I_k(1));
  I(2)=min(I(2),I_k(2));
  k=k+1;
end

if I(1)<I(2)       %%% Check that radii polynomials are really negative %%%
   r=(I(1)+I(2))/2;
   pr=p(:,4);
   for j=3:-1:1
     pr=p(:,j)+intval(r)*pr;
   end
   if max(sup(pr))>=0
      I=[1,-1];
   end
end

if I(1)>=I(2)    %%% empty interval %%%
    
    if dal0>0
        NEXTSTEP=1;  %%% Already had a successful step %%%
    else
        if (9/10)*dal >= dal_min
            dal=(9/10)*dal;        
        else
            %%% (Failure of the continuation step) %%%
            display(['Failure! Problem with finding dal at alpha = ',num2str(alpha0)])
            B0=B0_intermediate(:,count);
            return
        end
    end
    
else %%% I is not empty %%%
   
    %%% Store successful parameters %%%
    dal0=dal;
    r0=r;
    
    if (10/9)*dal<=dal_max       
        dal=(10/9)*dal;    
    else 
        NEXTSTEP=1;
    end
    
end  %%% End of the if Empty/Nonempty(I) %%%

end  %%% End of the while NEXTSTEP==0 %%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%                                    %%%
%%%             NEXT STEP              %%%
%%% (Success of the continuation step) %%%
%%%    We first save and then update   %%%
%%%                                    %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

count=count+1;

alpha0=inf(ialpha0+intval(dal0));
hx=bx+dal0*dx;
dal=dal0;

B0_intermediate(:,count+1)=[alpha0;r0;hx];

if max(abs(hx(2*m-1:2*m)))>1.e-9
    change_dimension=1;
end

display(['count='  num2str(count),'   dal0=' num2str(dal0),'  alpha0=' num2str(alpha0),'  r0=' num2str(r0),'  m=' num2str(m),'  pi/2-alpha0=' num2str(alpha0-pi/2)])

dal0=0;
    
end  %%% End of the "while(alpha0<alpha_max)" %%%

B0=B0_intermediate(:,count+1);
B31=B0_intermediate(:,count+1);
save B31 B31
disp('Success!');

end

function [F]=F_Wright(x,alpha)
m=length(x)/2;
L=x(1,1);
ba=[x(2,1);x((2*(1:m-1)'+1),1)];
bb=[0;x((2*(2:m)'),1)];

F(1,1)=ba(1,1)+2*sum(ba((2:m),1));
F(2,1)=alpha*ba(1,1)*(1+ba(1,1))+2*alpha*sum(cos(L*(1:m-1))'.*(ba((2:m),1).^2+bb((2:m),1).^2));

for k=1:m-1  
    R=[[alpha*cos(k*L) -k*L+alpha*sin(k*L)];[k*L-alpha*sin(k*L) alpha*cos(k*L)]];
    Nonlinear=[0;0];
    for k1=-m+1:m-1
        k2=k-k1;
        if abs(k2)<m
            Theta=[[cos(k1*L) sin(k1*L)];[-sin(k1*L) cos(k1*L)]];
            Nonlinear=Theta*[a(ba,k1)*a(ba,k2)-b(bb,k1)*b(bb,k2);a(ba,k1)*b(bb,k2)+b(bb,k1)*a(ba,k2)]+Nonlinear;
        end
    end
    F([2*k+1;2*k+2],1)=R*[a(ba,k);b(bb,k)]+alpha*Nonlinear;
end

end

function [J]=F_Wright_x(x,alpha)

m=length(x)/2;
L=x(1,1);
ba=[x(2,1);x((2*(1:m-1)'+1),1)];
bb=[0;x((2*(2:m)'),1)];

tempJ1=[0 1 2*ones(1,2*m-2)]; tempJ1(1,(2*(2:m)))=zeros(m-1,1);
J(1,:)=tempJ1*L/L;

J(2,1)=-2*alpha*sum((1:m-1)'.*sin(L*(1:m-1))'.*(ba((2:m),1).^2+bb((2:m),1).^2));
J(2,2)=alpha*(1+2*ba(1,1));

temp(1,(2*(1:m-1)))=(1:m-1); temp(1,(2*(0:m-2)+1))=(1:m-1);
J(2,(3:2*m))=4*alpha*cos(L*temp).*x((3:2*m),1)';

for k=1:m-1  

    R_L=[[-alpha*k*sin(k*L) -k+alpha*k*cos(k*L)];[k-alpha*k*cos(k*L) -alpha*k*sin(k*L)]];
    Nonlinear=[0;0];
    for k1=-m+1:m-1
        k2=k-k1;
        if abs(k2)<m
            Theta_L=[[-k1*sin(k1*L) k1*cos(k1*L)];[-k1*cos(k1*L) -k1*sin(k1*L)]];
            Nonlinear=Theta_L*[a(ba,k1)*a(ba,k2)-b(bb,k1)*b(bb,k2);a(ba,k1)*b(bb,k2)+b(bb,k1)*a(ba,k2)]+Nonlinear;
        end
    end
    J([2*k+1;2*k+2],1)=R_L*[a(ba,k);b(bb,k)]+alpha*Nonlinear;
    
    J([2*k+1;2*k+2],2)=alpha*(Theta(0,L)+Theta(k,L))*[a(ba,k);b(bb,k)];    
    
    R=[[alpha*cos(k*L) -k*L+alpha*sin(k*L)];[k*L-alpha*sin(k*L) alpha*cos(k*L)]];
    
    for l=1:m-1
        tempJ_kl=isequal(k,l)*R+alpha*(Theta(l,L)+Theta(k-l,L))*[[a(ba,k-l) -b(bb,k-l)];[b(bb,k-l) a(ba,k-l)]];
        if k+l>=m
            J([2*k+1;2*k+2],[2*l+1;2*l+2])=tempJ_kl;
        else
            J([2*k+1;2*k+2],[2*l+1;2*l+2])=tempJ_kl+alpha*(Theta(-l,L)+Theta(k+l,L))*[[a(ba,k+l) b(bb,k+l)];[b(bb,k+l) -a(ba,k+l)]];
        end
    end
    
end

end 

function [F]=F_Wright_alpha(x)
m=length(x)/2;
L=x(1,1);
ba=[x(2,1);x((2*(1:m-1)'+1),1)];
bb=[0;x((2*(2:m)'),1)];

F(1,1)=0*L;
F(2,1)=ba(1,1)*(1+ba(1,1))+2*sum(cos(L*(1:m-1))'.*(ba((2:m),1).^2+bb((2:m),1).^2));

for k=1:m-1  
    R_alpha=[[cos(k*L) sin(k*L)];[-sin(k*L) cos(k*L)]];
    Nonlinear=[0;0];
    for k1=-m+1:m-1
        k2=k-k1;
        if abs(k2)<m
            Theta=[[cos(k1*L) sin(k1*L)];[-sin(k1*L) cos(k1*L)]];
            Nonlinear=Theta*[a(ba,k1)*a(ba,k2)-b(bb,k1)*b(bb,k2);a(ba,k1)*b(bb,k2)+b(bb,k1)*a(ba,k2)]+Nonlinear;
        end
    end
    F([2*k+1;2*k+2],1)=R_alpha*[a(ba,k);b(bb,k)]+Nonlinear;
end

end

function [T]=Theta(k,L)
     T=[[cos(k*L) sin(k*L)];[-sin(k*L) cos(k*L)]];
end

function [coeff]=a(ba,n1)
    m=length(ba);
    if abs(n1)<m
        coeff=ba(abs(n1)+1,1);
    else
        coeff=0;
    end
end

function [coeff]=b(bb,n1)
    m=length(bb);
    if abs(n1)<m
        if n1<0
            coeff=-(bb(abs(n1)+1,1));
        else
            coeff=bb(n1+1,1);
        end
    else
        coeff=0;
    end
end

function [w]=omega(k,s)

if k==0
    w=1^s;
else 
    w=abs(k)^s;
end 

end

function [positive_roots]=rootsfinder_radii_polynomial(coeff)
% finds the positive roots of a radii polynomial;
% if no roots are found, it returns [1,-1]. 

tol=1e-3;
maxiter=10;
positive_roots=[1,-1];

D=coeff(2)^2-4*coeff(1)*coeff(3);
if D>0 && coeff(2)<0 
  r1=(-coeff(2)-sqrt(D))/(2*coeff(3));
  [root(1),froot(1)]=newtonpoly(coeff,r1,tol,maxiter);
  if root(1)>=0 && abs(froot(1))<tol
    r2=(-coeff(2)+sqrt(D))/(2*coeff(3));
    [root(2),froot(2)]=newtonpoly(coeff,r2,10*tol,maxiter);
    if root(2)>root(1) && abs(froot(2))<10*tol 
      positive_roots=root;   
    end
  end
end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [root,froot]=newtonpoly(coeff,x0,tol,maxiter)
% Newton iterations for a polynomial

n=length(coeff)-1;
coefdf=coeff(2:n+1).*(1:n);
root=x0;
itercount=0;
diff=abs(root)*tol+1;
while abs(diff)>abs(root)*tol && itercount<maxiter,
   froot=coeff(n+1);
   for j=n:-1:1
      froot=coeff(j)+root*froot;
   end
   dfroot=coefdf(n);
   for j=n-1:-1:1
      dfroot=coefdf(j)+root*dfroot;
   end
   if dfroot~=0
      diff=-froot/dfroot;
      root=root+diff;    
   else
     diff=0;
   end
   itercount=itercount+1;
end

end