%---------------------------------------------------|
%----MULTI REGIME - SECOND ORDER - FINITE BUFFER----|
%---------WITH ORDERED SCHUR DECOMPOSITION----------|
%-----------NONE OF THE RATES IS ZERO---------------|
%------  NONE OF THE STATES IS REPULSIVE------------|
%---------------------------------------------------|

%---------Output: Statewise PDF Plot----------------|
%---------Output: Total PDF Plot--------------------|
%---------Output: Statewise PDF Plot----------------|
%---------Output: Total CDF Plot--------------------|

clear;
%------------------PARAMETERS-----------------------|                                                   
%thresholds=[0,3,10];
%driftsR{1}=[4,3,-2,-1];%drifts in regime 1
%driftsR{2}=[4,3,-2,-1];%drifts in regime 2
%driftsT{1}=[4,3,-2,-1];%drifts in regime 1
%driftsT{2}=[4,3,-2,-1];%drifts in regime 1
%driftsT{3}=[4,3,-2,-1];%drifts in regime 1
%varR{1}=[1,1,1,1];
%varR{2}=[1,0,1,1];
%varT{1}=[1,1,1,1];
%varT{2}=[1,0,1,1];
%varT{3}=[1,0,1,1];
%RefvsAbs{1}=[1,1,1,1];%at T^(0)=0; 1 for reflecting, 0 for absorbing, 2 for first order 
%RefvsAbs{2}=[1,0,1,1];%at T^(K)
%QR{1}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];       %|
%QR{2}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];       %|
%QT{1}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];       %|
%QT{2}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];       %|
%QT{3}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];       %|
thresholds=[0,0.5,1];
driftsR{1}=[1,2];%drifts in regime 1
driftsR{2}=[-1,2];%drifts in regime 2
driftsT{1}=[1,2];%drifts in regime 1
driftsT{2}=[0,2];%drifts in regime 1
driftsT{3}=[-1,0];%drifts in regime 1

varR{1}=[0.1,0.1]*10;
varR{2}=[0.1,0.1]*10;
varT{1}=[0.1,0.1]*10;
varT{2}=[0,0.1]*10;
varT{3}=[0.1,0.1]*10;
RefvsAbs{1}=[1,1];%at T^(0)=0; 1 for reflecting, 0 for absorbing, 2 for first order 
RefvsAbs{2}=[1,1];%at T^(K)
QR{1}=[-1 1;1,-1];       %|
QR{2}=[-1 1;1,-1];       %|
QT{1}=[-1 1;1,-1];       %|
QT{2}=[-1 1;1,-1];       %|
QT{3}=[-1 1;1,-1];       %|
noofsamples=1000;%For Plotting                     %|
%---------------------------------------------------|


%-------------------PREPARE-------------------------|
r=length(driftsR);                                  %|
s=length(driftsR{1});                               %|
%---------------SYSTEM RECOGNIZATION----------------|
%Classify intermediate boundary points
it11no=zeros(1,r-1);%number of type-1 intermediate boundary points
it12no=zeros(1,r-1);%number of type-1 intermediate boundary points
it2no=zeros(1,r-1);
it3no=zeros(1,r-1);
it4no=zeros(1,r-1);
it5no=zeros(1,r-1);
it6no=zeros(1,r-1);
it7no=zeros(1,r-1);
it8no=zeros(1,r-1);
it9no=zeros(1,r-1);
it10no=zeros(1,r-1);
itmatrix=zeros(r-1,s);
Catino=zeros(1,r-1);
for i=1:r-1
    for j=1:s
        if varR{i}(j)>0%+
            if varR{i+1}(j)>0
                if driftsR{i}(j)>0&driftsR{i+1}(j)<0
                    it11no(i)=it11no(i)+1;
                    it11{i}{it11no(i)}=j;
                    itmatrix(i,j)=11;                
                    Catino(i)=Catino(i)+1;
                    nonzeroCati{i}{Catino(i)}=j;          
                else
                    it12no(i)=it12no(i)+1;
                    it12{i}{it12no(i)}=j;
                    itmatrix(i,j)=12;                
                end
            elseif driftsR{i+1}(j)<0&varR{i+1}(j)==0
                it2no(i)=it2no(i)+1;
                it2{i}{it2no(i)}=j;                
                itmatrix(i,j)=2;                
            elseif driftsR{i+1}(j)>0&varR{i+1}(j)==0
                it4no(i)=it4no(i)+1;
                it4{i}{it4no(i)}=j;                
                itmatrix(i,j)=4;
            else
                it10no(i)=it10no(i)+1;
                it10{i}{it10no(i)}=j;                
                itmatrix(i,j)=10;
                Catino(i)=Catino(i)+1;
                nonzeroCati{i}{Catino(i)}=j;
            end
        elseif driftsR{i}(j)>0&varR{i}(j)==0%0+
            if varR{i+1}(j)>0
                it5no(i)=it5no(i)+1;
                it5{i}{it5no(i)}=j;                
                itmatrix(i,j)=5;                
            elseif driftsR{i+1}(j)>0&varR{i+1}(j)==0
                it8no(i)=it8no(i)+1;
                it8{i}{it8no(i)}=j;                
                itmatrix(i,j)=8;                
            else
                it10no(i)=it10no(i)+1;                
                it10{i}{it10no(i)}=j;
                itmatrix(i,j)=10;
                Catino(i)=Catino(i)+1;
                nonzeroCati{i}{Catino(i)}=j;
            end
        elseif driftsR{i}(j)<0&varR{i}(j)==0%0-
            if (driftsT{i+1}(j)==0|driftsT{i+1}(j)>0)&varR{i+1}(j)==0
                it3no(i)=it3no(i)+1;
                it3{i}{it3no(i)}=j;                
                itmatrix(i,j)=3;
                Catino(i)=Catino(i)+1;
                nonzeroCati{i}{Catino(i)}=j;
            elseif varR{i+1}(j)>0
                it7no(i)=it7no(i)+1;                
                it7{i}{it7no(i)}=j;
                itmatrix(i,j)=7;
            elseif driftsR{i+1}(j)<0&varR{i+1}(j)==0
                it9no(i)=it9no(i)+1;
                it9{i}{it9no(i)}=j;                
                itmatrix(i,j)=9;
            elseif (driftsT{i+1}(j)==0|driftsT{i+1}(j)<0)&varR{i+1}(j)==0
                it6no(i)=it6no(i)+1;
                it6{i}{it6no(i)}=j;                
                itmatrix(i,j)=6;
                Catino(i)=Catino(i)+1;
                nonzeroCati{i}{Catino(i)}=j;
            else
                it10no(i)=it10no(i)+1;
                it10{i}{it10no(i)}=j;                
                itmatrix(i,j)=10;
                Catino(i)=Catino(i)+1;
                nonzeroCati{i}{Catino(i)}=j;
            end
        end
    end
end
%Classify terminal boundary points
t0t1no=0;%number of type-1 terminal 0 boundary points
t0t2no=0;
t0t3no=0;
t0tmatrix=zeros(1,s);
Cat0no=0;
for j=1:s
    if RefvsAbs{1}(j)==1%reflecting        
        t0t1no=t0t1no+1;
        t0t1{t0t1no}=j;
        t0tmatrix(j)=1;        
    elseif varR{1}(j)==0&driftsR{1}(j)>0
        t0t1no=t0t1no+1;
        t0t1{t0t1no}=j;
        t0tmatrix(j)=1;        
    elseif RefvsAbs{1}(j)==0&varR{1}(j)>0%absorbing second order
        t0t2no=t0t2no+1;
        t0t2{t0t2no}=j;
        t0tmatrix(j)=2;        
        Cat0no=Cat0no+1;
        nonzeroCat0{Cat0no}=j;
    else
        t0t3no=t0t3no+1;
        t0t3{t0t3no}=j;
        t0tmatrix(j)=3;
        Cat0no=Cat0no+1;
        nonzeroCat0{Cat0no}=j;
    end        
end

tTKt1no=0;%number of type-1 terminal T^{K} boundary points
tTKt2no=0;
tTKt3no=0;
tTKtmatrix=zeros(1,s);
CatTKno=0;
for j=1:s
    if RefvsAbs{2}(j)==1%reflecting        
        tTKt1no=tTKt1no+1;
        tTKt1{tTKt1no}=j;
        tTKtmatrix(j)=1;        
    elseif varR{r}(j)==0&driftsR{r}(j)<0
        tTKt1no=tTKt1no+1;
        tTKt1{tTKt1no}=j;
        tTKtmatrix(j)=1;        
    elseif RefvsAbs{2}(j)==0&varR{r}(j)>0%absorbing second order
        tTKt2no=tTKt2no+1;
        tTKt2{tTKt2no}=j;
        tTKtmatrix(j)=2;      
        CatTKno=CatTKno+1;
        nonzeroCatTK{CatTKno}=j;
    else
        tTKt3no=tTKt3no+1;
        tTKt3{tTKt3no}=j;
        tTKtmatrix(j)=3;
        CatTKno=CatTKno+1;
        nonzeroCatTK{CatTKno}=j;
    end        
end


%---------------------------------------------------|


%Reorder states:
%[f+ f0 f*]
reorder=zeros(r,s);%+,0,*
counts=ones(1,r);
countsplus=zeros(1,r);
countszero=zeros(1,r);
countsstar=zeros(1,r);
for i=1:r
    for j=1:s
        if varR{i}(j)>0%+
            newvarR{i}(counts(i),counts(i))=varR{i}(j);
            newdriftsR{i}(counts(i),counts(i))=driftsR{i}(j);
            reorder(i,j)=counts(i);
            counts(i)=counts(i)+1;            
            countsplus(i)=countsplus(i)+1;            
        end
    end
end
for i=1:r
    for j=1:s
        if varR{i}(j)==0&driftsR{i}(j)==0%0
            newvarR{i}(counts(i),counts(i))=varR{i}(j);
            newdriftsR{i}(counts(i),counts(i))=driftsR{i}(j);
            reorder(i,j)=counts(i);
            counts(i)=counts(i)+1;
            countszero(i)=countszero(i)+1;
        end
    end
end
for i=1:r
    for j=1:s
        if varR{i}(j)==0&(driftsR{i}(j)>0|driftsR{i}(j)<0)%*
            newvarR{i}(counts(i),counts(i))=varR{i}(j);
            newdriftsR{i}(counts(i),counts(i))=driftsR{i}(j);
            reorder(i,j)=counts(i);
            counts(i)=counts(i)+1;
            countsstar(i)=countsstar(i)+1;
        end
    end
end
%Construct newQ matrices
for i=1:r
    for j=1:s
        for k=1:s
            newQR{i}(j,k)=QR{i}(reorder(i,j),reorder(i,k));
        end
    end
end
for i=1:r
    Qplusplus{i}=newQR{i}(1:countsplus(i),1:countsplus(i));
    Qpluszero{i}=newQR{i}(1:countsplus(i),countsplus(i)+1:countsplus(i)+countszero(i));
    Qplusstar{i}=newQR{i}(1:countsplus(i),countsplus(i)+countszero(i)+1:s);
    Qzeroplus{i}=newQR{i}(countsplus(i)+1:countsplus(i)+countszero(i),1:countsplus(i));
    Qzerozero{i}=newQR{i}(countsplus(i)+1:countsplus(i)+countszero(i),countsplus(i)+1:countsplus(i)+countszero(i));
    Qzerostar{i}=newQR{i}(countsplus(i)+1:countsplus(i)+countszero(i),countsplus(i)+countszero(i)+1:s);
    Qstarplus{i}=newQR{i}(countsplus(i)+countszero(i)+1:s,1:countsplus(i));
    Qstarzero{i}=newQR{i}(countsplus(i)+countszero(i)+1:s,countsplus(i)+1:countsplus(i)+countszero(i));
    Qstarstar{i}=newQR{i}(countsplus(i)+countszero(i)+1:s,countsplus(i)+countszero(i)+1:s);
    varplus{i}=newvarR{i}(1:countsplus(i),1:countsplus(i));
    driftsplus{i}=newdriftsR{i}(1:countsplus(i),1:countsplus(i));    
    driftsstar{i}=newdriftsR{i}(countsplus(i)+1:s,countsplus(i)+1:s);
end

countsminus=zeros(1,r);
for i=1:r
    for j=1:s
        if driftsR{i}(j)<0&varR{i}(j)==0
            countsminus(i)=countsminus(i)+1;
        end
    end
end

%B)Ordered Schur Decomposition:                    %|
%------------------------------                    %|
for i=1:r       
    [U,U_inv,A,Aplus,Aminus,Nminus,pik,lambda] = additive_decompose_mixed(varplus{i}./2,driftsplus{i},driftsstar{i},Qplusplus{i},Qplusstar{i},Qstarplus{i},Qstarstar{i},countsminus(i));
    
    Nminusc{i}=Nminus;                             %|
    Nplusc{i} = 2*countsplus(i)+countsstar(i)-Nminus-1;                      %|
    Uinvc{i}=U_inv;                                %|
    Lzero = U_inv(Nminus+1,1:countsplus(i));                   %|
    Lminus = U_inv(1:Nminus,1:countsplus(i));                  %|
    Lplus = U_inv(Nminus+2:2*countsplus(i)+countsstar(i),1:countsplus(i));               %|
    Mzero = U_inv(Nminus+1,countsplus(i)+1:2*countsplus(i));               %|
    Mminus = U_inv(1:Nminus,countsplus(i)+1:2*countsplus(i));              %|
    Mplus = U_inv(Nminus+2:2*countsplus(i)+countsstar(i),countsplus(i)+1:2*countsplus(i));           %|
    Pzero = U_inv(Nminus+1,2*countsplus(i)+1:2*countsplus(i)+countsstar(i));               %|
    Pminus = U_inv(1:Nminus,2*countsplus(i)+1:2*countsplus(i)+countsstar(i));              %|
    Pplus = U_inv(Nminus+2:2*countsplus(i)+countsstar(i),2*countsplus(i)+1:2*countsplus(i)+countsstar(i));           %|    
    gecicizero=[Lzero Pzero];
    geciciminus=[Lminus Pminus];
    geciciplus=[Lplus Pplus];
    Lnewzero=zeros(1,s);
    Lnewminus=zeros(Nminus,s);
    Lnewplus=zeros(Nplusc{i},s);
    for j=1:s
        Lnewzero(j)=gecicizero(reorder(i,j));
        Lnewminus(:,j)=geciciminus(:,reorder(i,j));
        Lnewplus(:,j)=geciciplus(:,reorder(i,j));
    end
    gecicizero=[Mzero zeros(1,countsstar(i))];
    geciciminus=[Mminus Aminus*Pminus];
    geciciplus=[Mplus Aplus*Pplus];
    Mnewzero=zeros(1,s);
    Mnewminus=zeros(Nminus,s);
    Mnewplus=zeros(Nplusc{i},s);
    for j=1:s
        Mnewzero(j)=gecicizero(reorder(i,j));
        Mnewminus(:,j)=geciciminus(:,reorder(i,j));
        Mnewplus(:,j)=geciciplus(:,reorder(i,j));
    end
        
    
    Lzeroc{i}=Lnewzero;                               %|
    Lplusc{i}=Lnewplus;                               %|
    Lminusc{i}=Lnewminus;                             %|
    Mzeroc{i}=Mnewzero;                               %|
    Mplusc{i}=Mnewplus;                               %|
    Mminusc{i}=Mnewminus;                             %|
    
    
    
    
    
    Aplusc{i}=Aplus;                               %|
    Aminusc{i}=Aminus;                             %|  
end                                                %|
%---------------------------------------------------|

%---------------------------------------------------|



%---------------------SOLUTION----------------------|
rownos=zeros(1,r+1);
for i=2:r+1
    rownos(i)=2*sum(countsplus(1:i-1))+sum(countsstar(1:i-1));
end
rows=2*sum(countsplus)+sum(countsstar)+Cat0no+CatTKno+sum(Catino);
columns=(r-1)*s+sum(it12no)+sum(it3no)+sum(it4no)+sum(it6no)+sum(it7no)+2*sum(it11no)+s+t0t2no+s+tTKt2no+1;

totalmatrix=zeros(rows,columns);%                  %|



%A) Boundary Conditions:                           %|
%-----------------------                           %|
%Boundary C. at Zero                               %| 
y1=diag(driftsR{1});                                %|
z1=diag(varR{1}./2);                          %|

%T-1
nextstep1=1;%starting column no
nextstep2=s;%ending column no
totalmatrix(rownos(1)+1:rownos(2),nextstep1:nextstep2)=[Lzeroc{1}*y1-Mzeroc{1}*z1;...
 Lminusc{1}*y1-Mminusc{1}*z1;...                   %|
 expm(-Aplusc{1}*thresholds(2))*(Lplusc{1}*y1-Mplusc{1}*z1)]; 

%T-5
if t0t2no>0
    nextstep1=nextstep2+1;
    nextstep2=nextstep2+t0t2no;
    gecici=zeros(rownos(2),t0t2no);
    gecici2=[Lzeroc{1};Lminusc{1};expm(-Aplusc{1}*thresholds(2))*Lplusc{1}];
    for i=1:t0t2no
        gecici(:,i)=gecici2(:,t0t2{i});        
    end    
    totalmatrix(rownos(1)+1:rownos(2),nextstep1:nextstep2)=gecici;
end
%Boundary C. at Intermediate boundary points                         %| 
%I-1
for i=1:r-1
    y1=diag(driftsR{i});                        %|
    z1=diag(varR{i}./2);                  %|
    y2=diag(driftsR{i+1});                      %|
    z2=diag(varR{i+1}./2);                %|
    nextstep1=nextstep2+1;
    nextstep2=nextstep2+s;    
    totalmatrix(rownos(i)+1:rownos(i+1),nextstep1:nextstep2)=[-Lzeroc{i}*y1+Mzeroc{i}*z1;...
         -expm(Aminusc{i}*(thresholds(i+1)-thresholds(i)))*(Lminusc{i}*y1-Mminusc{i}*z1);...
         -Lplusc{i}*y1+Mplusc{i}*z1];              %|
    totalmatrix(rownos(i+1)+1:rownos(i+2),nextstep1:nextstep2)=[Lzeroc{i+1}*y2-Mzeroc{i+1}*z2;...
         Lminusc{i+1}*y2-Mminusc{i+1}*z2;...       %|
         expm(-Aplusc{i+1}*(thresholds(i+2)-thresholds(i+1)))*(Lplusc{i+1}*y2-Mplusc{i+1}*z2)];
     
end
%I-2
for i=1:r-1
    if it12no(i)>0
        nextstep1=nextstep2+1;
        nextstep2=nextstep2+it12no(i);
        z1=diag(varR{i}./2);                  %|        
        z2=diag(varR{i+1}./2);                %|       
        gecici2=[Lzeroc{i}*z1;...
         +expm(Aminusc{i}*(thresholds(i+1)-thresholds(i)))*Lminusc{i}*z1;...
         +Lplusc{i}*z1];              %|
        gecici3=zeros(rownos(i+1)-rownos(i),it12no(i));
        for j=1:it12no(i)
            gecici3(:,j)=gecici2(:,it12{i}{j});
        end    
        totalmatrix(rownos(i)+1:rownos(i+1),nextstep1:nextstep2)=gecici3;                
        gecici2=[-Lzeroc{i+1}*z2;...
         -Lminusc{i+1}*z2;...       %|
         -expm(-Aplusc{i+1}*(thresholds(i+2)-thresholds(i+1)))*Lplusc{i+1}*z2];     
        gecici3=zeros(rownos(i+2)-rownos(i+1),it12no(i));
        for j=1:it12no(i)
            gecici3(:,j)=gecici2(:,it12{i}{j});
        end            
        totalmatrix(rownos(i+1)+1:rownos(i+2),nextstep1:nextstep2)=gecici3;
    end
end
%I-4
for i=1:r-1
    if it3no(i)+it4no(i)+it11no(i)>0
        nextstep1=nextstep2+1;
        nextstep2=nextstep2+it3no(i)+it4no(i)+it11no(i);                             
        gecici2=[Lzeroc{i};...
         +expm(Aminusc{i}*(thresholds(i+1)-thresholds(i)))*Lminusc{i};...
         +Lplusc{i}];              %|
        gecici3=zeros(rownos(i+1)-rownos(i),it3no(i)+it4no(i)+it11no(i));
        gecicicount=1;
        for j=1:it3no(i)
            gecici3(:,j)=gecici2(:,it3{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                    
        gecicicount=1;
        for j=it3no(i)+1:it3no(i)+it4no(i)
            gecici3(:,j)=gecici2(:,it4{i}{gecicicount});
            gecicicount=gecicicount+1;
        end           
        gecicicount=1;
        for j=it3no(i)+it4no(i)+1:it3no(i)+it4no(i)+it11no(i)
            gecici3(:,j)=gecici2(:,it11{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                            
        
        totalmatrix(rownos(i)+1:rownos(i+1),nextstep1:nextstep2)=gecici3;
    end
end
%I-6
for i=1:r-1
    if it6no(i)+it7no(i)+it11no(i)>0
        nextstep1=nextstep2+1;
        nextstep2=nextstep2+it6no(i)+it7no(i)+it11no(i);                                     
        gecici2=[Lzeroc{i+1};...
         Lminusc{i+1};...       %|
         expm(-Aplusc{i+1}*(thresholds(i+2)-thresholds(i+1)))*Lplusc{i+1}];     
        gecici3=zeros(rownos(i+2)-rownos(i+1),it6no(i)+it7no(i)+it11no(i));
        gecicicount=1;
        for j=1:it6no(i)
            gecici3(:,j)=gecici2(:,it6{i}{gecicicount});
            gecicicount=gecicicount+1;
        end           
        gecicicount=1;
        for j=it6no(i)+1:it6no(i)+it7no(i)
            gecici3(:,j)=gecici2(:,it7{i}{gecicicount});
            gecicicount=gecicicount+1;
        end           
        gecicicount=1;
        for j=it6no(i)+it7no(i)+1:it6no(i)+it7no(i)+it11no(i)
            gecici3(:,j)=gecici2(:,it11{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                            
        
        totalmatrix(rownos(i+1)+1:rownos(i+2),nextstep1:nextstep2)=gecici3;
    end
end

%T-3
nextstep1=nextstep2+1;
nextstep2=nextstep2+s;                      
y1=diag(driftsR{r});                        %|
z1=diag(varR{r}./2);                  %|
totalmatrix(rownos(r)+1:rownos(r+1),nextstep1:nextstep2)=[-Lzeroc{r}*y1+Mzeroc{r}*z1;...
         -expm(Aminusc{r}*(thresholds(r+1)-thresholds(r)))*(Lminusc{r}*y1-Mminusc{r}*z1);...
         -Lplusc{r}*y1+Mplusc{r}*z1];              %|
%T-6
if tTKt2no>0
    nextstep1=nextstep2+1;
    nextstep2=nextstep2+tTKt2no;
    gecici2=[Lzeroc{r};...
         +expm(Aminusc{r}*(thresholds(r+1)-thresholds(r)))*Lminusc{r};...
         +Lplusc{r}];              %|
    
    gecici3=zeros(2*countsplus(r)+countsstar(r),tTKt2no);
    for i=1:tTKt2no
        gecici3(:,i)=gecici2(:,tTKt2{i});
    end    
    totalmatrix(rownos(r)+1:rownos(r+1),nextstep1:nextstep2)=gecici3;
end

%add terms containing C
%T1-I1-T3

%T1
nextstepc1=1;%starting column no
nextstepc2=s;%ending column no
nextstepr1=2*sum(countsplus)+sum(countsstar);
nextstepr2=2*sum(countsplus)+sum(countsstar);
if Cat0no>0
    nextstepr1=2*sum(countsplus)+sum(countsstar)+1;%starting row no
    nextstepr2=2*sum(countsplus)+sum(countsstar)+Cat0no;%ending row no
    gecici=zeros(nextstepr2-nextstepr1+1,nextstepc2-nextstepc1+1);
    for i=1:nextstepr2-nextstepr1+1
        gecici(i,:)=-QT{1}(nonzeroCat0{i},:);
    end
    totalmatrix(nextstepr1:nextstepr2,nextstepc1:nextstepc2)=gecici;
end

%T5
if t0t2no>0
    nextstepc1=nextstepc2+1;
    nextstepc2=nextstepc2+t0t2no;
end

%I1
for i=1:r-1
    nextstepc1=nextstepc2+1;
    nextstepc2=nextstepc2+s;
    if Catino(i)>0
        nextstepr1=nextstepr2+1;%starting row no
        nextstepr2=nextstepr2+Catino(i);%ending row no
        gecici=zeros(nextstepr2-nextstepr1+1,nextstepc2-nextstepc1+1);
        for j=1:nextstepr2-nextstepr1+1
            gecici(j,:)=-QT{i+1}(nonzeroCati{i}{j},:);
        end
        totalmatrix(nextstepr1:nextstepr2,nextstepc1:nextstepc2)=gecici;
    end
end
%I2
for i=1:r-1
    if it12no(i)>0
        nextstepc1=nextstepc2+1;
        nextstepc2=nextstepc2+it12no(i);
    end
end
%I4
for i=1:r-1
    if it3no(i)+it4no(i)+it11no(i)>0
        nextstepc1=nextstepc2+1;
        nextstepc2=nextstepc2+it3no(i)+it4no(i)+it11no(i);                             
    end
end
%I6
for i=1:r-1
    if it6no(i)+it7no(i)+it4no(i)>0
        nextstepc1=nextstepc2+1;
        nextstepc2=nextstepc2+it6no(i)+it7no(i)+it11no(i);                             
    end
end
%T3
nextstepc1=nextstepc2+1;
nextstepc2=nextstepc2+s;
if CatTKno>0    
    nextstepr1=nextstepr2+1;%starting row no
    nextstepr2=nextstepr2+CatTKno;%ending row no
    gecici=zeros(nextstepr2-nextstepr1+1,nextstepc2-nextstepc1+1);
    for i=1:nextstepr2-nextstepr1+1
        gecici(i,:)=-QT{r+1}(nonzeroCatTK{i},:);
    end
    totalmatrix(nextstepr1:nextstepr2,nextstepc1:nextstepc2)=gecici;
end
%Normalization conditions
%for C unknowns
totalmatrix(2*sum(countsplus)+sum(countsstar)+1:2*sum(countsplus)+sum(countsstar)+Cat0no+CatTKno+sum(Catino),columns)=ones(Cat0no+CatTKno+sum(Catino),1);

%for regimes
for i=1:r
    integ=[Lzeroc{i}*(thresholds(i+1)-thresholds(i));...
         Aminusc{i}^-1*(expm(Aminusc{i}*(thresholds(i+1)-thresholds(i)))-eye(Nminusc{i}))*Lminusc{i};...
         Aplusc{i}^-1*(eye(Nplusc{i})-expm(-Aplusc{i}*(thresholds(i+1)-thresholds(i))))*Lplusc{i}]*ones(s,1);
    totalmatrix(rownos(i)+1:rownos(i+1),columns)=integ;
end       
%-   -   -   -   -   -   -   -   -   -   -   -   - %|





                    %|
%Boundary C. at Buffer Length                      %| 
%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%B)Solution:                                       %|
%-----------                                       %|
totalvector=[zeros(1,columns-1) 1];                    %|
solution=totalvector/totalmatrix;                  %|
%---------------------------------------------------|


%---------------------PLOTS-------------------------|
samples=2000;                                    %|
deltax=thresholds(r+1)/samples;                       %|
pdfs=zeros(s,samples);                             %|
cdfs=zeros(s,samples);                             %|
totalpdfs=zeros(1,samples);                        %|
totalcdfs=zeros(1,samples);                        %|


for i=1:r                                          %|
    azeroc{i}=solution(rownos(i)+1);
    aminusc{i}=solution(rownos(i)+2:rownos(i)+1+Nminusc{i});
    aplusc{i}=solution(rownos(i)+2+Nminusc{i}:rownos(i)+1+Nminusc{i}+Nplusc{i});
end                                                %|
i=1;                                               %|
for j=1:samples                                    %|
    if j*deltax<=thresholds(i+1)                   %|
        pdfs(:,j)=(azeroc{i}*Lzeroc{i}+...         %|
         aminusc{i}*expm(Aminusc{i}*((j-1)*deltax-thresholds(i)))*Lminusc{i}+...
         aplusc{i}*expm(-Aplusc{i}*(thresholds(i+1)-(j-1)*deltax))*Lplusc{i})';
    elseif i<r                                     %|
        i=i+1;                                     %|
        pdfs(:,j)=(azeroc{i}*Lzeroc{i}+...         %|
         aminusc{i}*expm(Aminusc{i}*((j-1)*deltax-thresholds(i)))*Lminusc{i}+...
         aplusc{i}*expm(-Aplusc{i}*(thresholds(i+1)-(j-1)*deltax))*Lplusc{i})';
    else                                           %|
        pdfs(:,j)=(azeroc{i}*Lzeroc{i}+...         %|
         aminusc{i}*expm(Aminusc{i}*((j-1)*deltax-thresholds(i)))*Lminusc{i}+...
         aplusc{i}*Lplusc{i})';                    %|
    end                                            %|
end                                                %|
t=0:deltax:deltax*(samples-1);                     %|
figure(1);                                         %|
hold on                                            %|
for i=1:s                                          %|
    plot(t,pdfs(i,:));                             %|
end                                                %|




i=1;                                               %|
m=0;                                               %|
for j=1:samples                                    %|
    if j*deltax<=thresholds(i+1)                   %|
        cdfs(:,j)= m+(azeroc{i}*Lzeroc{i}*((j-1)*deltax-thresholds(i))+...
         aminusc{i}*(Aminusc{i}^-1)*(expm(Aminusc{i}*((j-1)*deltax-thresholds(i)))-eye(Nminusc{i}))*Lminusc{i}+...
         aplusc{i}*Aplusc{i}^-1*(expm(-Aplusc{i}*(thresholds(i+1)-(j-1)*deltax))-expm(-Aplusc{i}*(thresholds(i+1)-thresholds(i))))*Lplusc{i})';  
    elseif i<r                                     %|
        m=cdfs(:,j-1);                             %|
        i=i+1;                                     %|
        cdfs(:,j)= m+(azeroc{i}*Lzeroc{i}*((j-1)*deltax-thresholds(i))+...
         aminusc{i}*(Aminusc{i}^-1)*(expm(Aminusc{i}*((j-1)*deltax-thresholds(i)))-eye(Nminusc{i}))*Lminusc{i}+...
         aplusc{i}*Aplusc{i}^-1*(expm(-Aplusc{i}*(thresholds(i+1)-(j-1)*deltax))-expm(-Aplusc{i}*(thresholds(i+1)-thresholds(i))))*Lplusc{i})';
    else                                           %|
        cdfs(:,j)= m+(azeroc{i}*Lzeroc{i}*((j-1)*deltax-thresholds(i))+...
         aminusc{i}*(Aminusc{i}^-1)*(expm(Aminusc{i}*((j-1)*deltax-thresholds(i)))-eye(Nminusc{i}))*Lminusc{i}+...
         aplusc{i}*Aplusc{i}^-1*(expm(-Aplusc{i}*(thresholds(i+1)-(j-1)*deltax))-expm(-Aplusc{i}*(thresholds(i+1)-thresholds(i))))*Lplusc{i})';    
    end                                            %|
end                                                %|   




%add accumulations

if Cat0no+CatTKno+sum(Catino)>0
    count=1;
    accs=solution(rownos(r+1)+1:columns-1);
    if Cat0no>0
        for i=1:Cat0no
            cdfs(nonzeroCat0{i},:)=cdfs(nonzeroCat0{i},:)+accs(count);
            count=count+1;
        end
    end
    for i=1:r-1
        if Catino(i)>0
            for j=1:Catino(i)
                cdfs(nonzeroCati{i}{j},floor(thresholds(i+1)/deltax)+1:samples)=cdfs(nonzeroCati{i}{j},floor(thresholds(i+1)/deltax)+1:samples)+accs(count);
                count=count+1;
            end
        end
    end
    if CatTKno>0
        for i=1:CatTKno
            cdfs(nonzeroCatTK{i},samples)=cdfs(nonzeroCatTK{i},samples)+accs(count);
            count=count+1;
        end
    end
        
end
t=0:deltax:deltax*(samples-1);                     %|
figure(1);                                         %|
hold on                                            %|
for i=1:s                                          %|
    plot(t,pdfs(i,:));                             %|
end                                                %|
t=0:deltax:deltax*(samples-1);                     %|
figure(2);                                         %|
hold on                                            %|
for i=1:s                                          %|
    plot(t,cdfs(i,:));                             %|
end                                                %|
totalcdfs=ones(1,s)*cdfs;                          %|
totalpdfs=ones(1,s)*pdfs;                          %|
t=0:deltax:deltax*(samples-1);                     %|
figure(3);                                         %|
hold on                                            %|
plot(t,totalpdfs);                                 %|
t=0:deltax:deltax*(samples-1);                     %|
hold on                                            %|
figure(4);plot(t,totalcdfs);                       %|
%---------------------------------------------------|