function x = NTAlgorithm2(A,B,q,X0,Y0,beta1,beta2,epsilon)

[m1,m2]=size(X0);
mu0=trace(X0*Y0)/m1;
r0=A*svec(X0)+B*svec(Y0)-q;
X=X0;
Y=Y0;
tau0=mu0;
tau=tau0;

ntilde=m1*(m1+1)/2;

tolerance=max([mu0*m1,norm(r0,2)]);

counter=0;

muk_1=trace(X0*Y0);

ratio=zeros(100,1);

while (tolerance>epsilon)
    rhs=-[tau*r0/tau0;svec(X)];
    W=NT(X,Y);
%    cond(W);
    v=SolLE(A,B,W,rhs);
    [m,n]=size(v);
    deltaXp=smat(v(1:m/2));
    deltaYp=smat(v(m/2+1:m));
    delta=norm(Symmetrization(W,deltaXp*deltaYp),'fro')/tau;
    denominator=sqrt(1+4*delta/(beta2-beta1))+1;
    alpha=2/denominator;
    Xhat=X+alpha*deltaXp;
    Yhat=Y+alpha*deltaYp;
    sigma=1-alpha;
    What=NT(Xhat,Yhat);
%    cond(What)
    Yhatinverse=Yhat\eye(m1);
%    Whatinverse=What\eye(m1);
%    secondorderterm=Kronecker(Whatinverse,Yhat)\(Kronecker(Whatinverse,deltaYp)*svec(deltaXp));
    secondorderterm=zeros(ntilde,n);
    rhs=[zeros(ntilde,1);sigma*tau*svec(Yhatinverse)-svec(Xhat)-secondorderterm];
%    size(rhs)
    v=SolLE(A,B,What,rhs);
    [m,n]=size(v);
    deltaXc=smat(v(1:m/2));
    deltaYc=smat(v(m/2+1:m));
    X=Xhat+deltaXc;
    Y=Yhat+deltaYc;
    tau=sigma*tau;
    r=A*svec(X)+B*svec(Y)-q;
    tolerance=max([trace(X*Y),norm(r,2)]);
    counter=counter+1;
    muk=trace(X*Y);
    ratio(counter)=muk/muk_1;
    muk_1=muk;
end

ratio1=ratio(1:counter,1);
x=[counter;ratio1;svec(X);svec(Y)];