function [U_inv,Aplus,Aminus,Nminus] = additive2_decompose_schur(R,Q);
% N is the size of all the matrices
N = length(Q);
e = ones(N,1);
I = eye(N,N);
Z = zeros(N,N);
%CALCULATE LAMBDA(MEAN DRIFT)
Q1=[Q ones(N,1)];
pik=[zeros(1,N) 1]/Q1;
lambda = pik*R*e;
%CHECK FOR LAMBDA==0
if(abs(lambda)<10^-5)
    disp('WARNING!!!-lambda is zero-WARNING!!!')
else
    %CALCULATE nMÝNUS&nPLUS
    Nminus=-1;
    Nplus=0;
    for i=1:N
        if R(i,i)>0
            Nminus=Nminus+1;
        else
            Nplus=Nplus+1;
        end
    end
    if lambda<0
        Nminus = Nminus+1;
        Nplus = Nplus-1;
    end
    %DECOMPOSITION
    A=Q*R^-1;
    [Q2,R2] = schur(A);
    Y2=zeros(1,N);
    [Q1,R1] =ordschur(Q2,R2,'lhp');
    [YY,II]=min(abs(diag(R1)));
    zeroplace=II;
    for i=1:N
        Y2(i)=(N-i+1)*2;
        if i<N&i>Nminus
            if ~(R1(i+1,i)==0)
                Y2(i)=Y2(i)-2;
            end
        elseif i<=Nminus
            if ~(R1(i+1,i)==0)
                Y2(i)=Y2(i)+2;
            end
        end
    end
    Y2(zeroplace)=(N-Nminus)*2+1;
    if R1(Nminus+1,Nminus+1)<0
        Y2(zeroplace)=Y2(zeroplace)-2;
    end
    [Q1,R1] =ordschur(Q1,R1,Y2);
    %DEFINE A11 ETC.
    A11 = R1(1:Nminus+1,1:Nminus+1);
    A12 = R1(1:Nminus+1,Nminus+2:N);
    A22 = R1(Nminus+2:N,Nminus+2:N);
    
    X1 = lyap(A11,-A22,A12);
    X11 = X1(1:Nminus,:);
    X12 = X1(Nminus+1,:);
    Aminus = A11(1:Nminus,1:Nminus);
    Aminuszero = A11(1:Nminus,Nminus+1);
    X2 = -Aminus\Aminuszero;
    temp = [eye(Nminus) X2 X11;zeros(1,Nminus) 1 X12;zeros(Nplus,Nminus+1) eye(Nplus)];
    temp_inv = [eye(Nminus) -X2 -X11+X2*X12;zeros(1,Nminus) 1 -X12;zeros(Nplus,Nminus+1) eye(Nplus)];
    U = Q1*temp;
    U_inv = temp_inv*Q1';
    Aplus = A22;
end