%---------------------------------------------------|
%----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,1];
driftsR{1}=[4,3,-3,-1];%drifts in regime 1
driftsR{2}=[0,3,-3,-1];%drifts in regime 2
driftsT{1}=[0,3,0,0];%drifts in regime 1
driftsT{2}=[0,3,-3,-1];%drifts in regime 1

varR{1}=[1,0,0,1];
varR{2}=[0,0,0,0];
varT{1}=[0,0,0,0];
varT{2}=[0,0,0,0];

RefvsAbs{1}=[0,0,0,0];%at T^(0)=0; 1 for reflecting, 0 for absorbing, 2 for first order 

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];       %|


drawmax=20;%For Plotting                     %|
%---------------------------------------------------|


%-------------------PREPARE-------------------------|
r=length(driftsR);                                  %|
s=length(driftsR{1});                               %|
%Reorder states in each regime
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
%Construct B matrices
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);
    
    Mplusstar{i}=Qplusstar{i}-Qpluszero{i}*(Qzerozero{i})^-1*Qzerostar{i};
    Mplusplus{i}=Qplusplus{i}-Qpluszero{i}*(Qzerozero{i})^-1*Qzeroplus{i};
    Mstarstar{i}=Qstarstar{i}-Qstarzero{i}*(Qzerozero{i})^-1*Qzerostar{i};
    Mstarplus{i}=Qstarplus{i}-Qstarzero{i}*(Qzerozero{i})^-1*Qzeroplus{i};
    B{i}=[zeros(countsplus(i),countsplus(i)) Mplusstar{i}*(newdriftsR{i}(countsplus(i)+countszero(i)+1:s,countsplus(i)+countszero(i)+1:s))^-1 -Mplusplus{i}*(newvarR{i}(1:countsplus(i),1:countsplus(i))./2)^-1;...
          zeros(countsstar(i),countsplus(i)) Mstarstar{i}*(newdriftsR{i}(countsplus(i)+countszero(i)+1:s,countsplus(i)+countszero(i)+1:s))^-1 -Mstarplus{i}*(newvarR{i}(1:countsplus(i),1:countsplus(i))./2)^-1;...
          eye(countsplus(i)) zeros(countsplus(i),countsstar(i)) newdriftsR{i}(1:countsplus(i),1:countsplus(i))*(newvarR{i}(1:countsplus(i),1:countsplus(i))./2)^-1];      
end

for i=1:r                                          %|
    %-   -   -   -   -   -   -   -   -   -   -   - %|
    %A)Find Eigenvalues and Eigenvectors:          %|
    %------------------------------------          %|
    [v,d]=eig(B{i}');                          %|
    vnew1{i}=zeros(2*countsplus(i)+countsstar(i),countsplus(i)+countsstar(i)+countszero(i));
    for j=1:2*countsplus(i)+countsstar(i)
        fiplus=v(1:countsplus(i),j)';
        fistar=v(countsplus(i)+1:countsplus(i)+countsstar(i),j)';
        fizero=-fiplus*Qpluszero{i}*(Qzerozero{i})^-1-fistar*Qstarzero{i}*(Qzerozero{i})^-1;       
        vnew1{i}(j,:)=[fiplus fizero fistar];        
    end
    dnew{i}=d;
    %reorder columns of eigenvectors to return to the original order of states
    vnew{i}=vnew1{i};
    for j=1:s
        vnew{i}(:,j)=vnew1{i}(:,reorder(i,j));
    end
    
    negeigvalcount=0;
    if i==r
        for j=1:2*countsplus(r)+countsstar(r)
            if dnew{r}(j,j)<-10^-10
                negeigvalcount=negeigvalcount+1;
                negeigval{negeigvalcount}=j;
            end
        end
        dnewgecici=zeros(negeigvalcount);
        vnewgecici=zeros(negeigvalcount,s);
        for j=1:negeigvalcount
            dnewgecici(j,j)=dnew{r}(negeigval{j},negeigval{j});
            vnewgecici(j,:)=vnew{r}(negeigval{j},:);
        end
        vnew{r}=vnewgecici;
        dnew{r}=dnewgecici;
    end
        
        
end                                                %|


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

%---------------SYSTEM RECOGNIZATION----------------|
%Classify intermediate boundary points
it1no=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);
it11no=zeros(1,r-1);
it12no=zeros(1,r-1);
it13no=zeros(1,r-1);
it31no=zeros(1,r-1);
it32no=zeros(1,r-1);
it61no=zeros(1,r-1);
it62no=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
                it1no(i)=it1no(i)+1;
                it1{i}{it1no(i)}=j;
                itmatrix(i,j)=1;                
            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% driftsR{i+1}(j)==0&varR{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;
            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 driftsR{i+1}(j)==0&varR{i+1}(j)==0
                it32no(i)=it32no(i)+1;
                it32{i}{it32no(i)}=j;                
                itmatrix(i,j)=32;
                Catino(i)=Catino(i)+1;
                nonzeroCati{i}{Catino(i)}=j;
            elseif driftsR{i+1}(j)>0&varR{i+1}(j)==0&driftsT{i+1}(j)>0
                it31no(i)=it31no(i)+1;
                it31{i}{it31no(i)}=j;                
                itmatrix(i,j)=31;
            elseif driftsR{i+1}(j)>0&varR{i+1}(j)==0&driftsT{i+1}(j)<0
                it61no(i)=it61no(i)+1;
                it61{i}{it61no(i)}=j;                
                itmatrix(i,j)=61;                
            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&driftsR{i+1}(j)>0&varR{i+1}(j)==0
                it13no(i)=it13no(i)+1;
                it13{i}{it13no(i)}=j;                
                itmatrix(i,j)=13;
                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
        else% driftsR{i}(j)==0&varR{i}(j)==0%0-
            if varR{i+1}(j)>0
                it12no(i)=it12no(i)+1;                
                it12{i}{it12no(i)}=j;
                itmatrix(i,j)=12;
                Catino(i)=Catino(i)+1;
                nonzeroCati{i}{Catino(i)}=j;
            elseif varR{i+1}(j)==0&driftsR{i+1}(j)>0
                it62no(i)=it62no(i)+1;                
                it62{i}{it62no(i)}=j;
                itmatrix(i,j)=62;
                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
%---------------------------------------------------|

rownos=zeros(1,r+1);

for i=2:r
    rownos(i)=2*sum(countsplus(1:i-1))+sum(countsstar(1:i-1));
end
rownos(r+1)=rownos(r)+negeigvalcount;


%---------------------SOLUTION----------------------|


rows=2*sum(countsplus(1:r-1))+sum(countsstar(1:r-1))+Cat0no+sum(Catino)+negeigvalcount;

%I1->(r-1)*s
%I2->sum(it1no)
%I3->sum(it2no) iptal
%I4->sum(it3no)+sum(it4no)+sum(it11no)%sum(it4no)+sum(it11no)+sum(it13no)+s
%um(it31no)+sum(it32no)
%I5->sum(it5no) iptal
%I6->sum(it6no)+sum(it7no)+sum(it12no)%sum(it7no)+sum(it12no)+sum(it13no)su
%m(it61no)+sum(it62no)
%T1->s
%T5->t0t2no
%T3->s
%T6->tTKt2no
%1 for normalization condition
columns=(r-1)*s+sum(it1no)+sum(it4no)+sum(it11no)+sum(it13no)+sum(it31no)+sum(it32no)+sum(it7no)+sum(it12no)+sum(it13no)+sum(it61no)+sum(it62no)+s+t0t2no+1;
totalmatrix=zeros(rows,columns);%                  %|

%A) Boundary Conditions:                           %|
%-----------------------                           %|
%Boundary C. at Zero                               %| 
%ASSUME R>1
%T-1
nextstep1=1;%starting column no
nextstep2=s;%ending column no
totalmatrix(rownos(1)+1:rownos(2),nextstep1:nextstep2)=-dnew{1}*vnew{1}*diag(varR{1})./2+vnew{1}*diag(driftsR{1});

%T-5
if t0t2no>0
    nextstep1=nextstep2+1;
    nextstep2=nextstep2+t0t2no;
    gecici=zeros(rownos(2),t0t2no)
    for i=1:t0t2no
        gecici(:,i)=vnew{1}(:,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
    nextstep1=nextstep2+1;
    nextstep2=nextstep2+s;
    gecici=zeros(rownos(i+1)-rownos(i));
    for j=1:rownos(i+1)-rownos(i)
        gecici(j,j)=exp(dnew{i}(j,j)*thresholds(i+1));
    end
    totalmatrix(rownos(i)+1:rownos(i+1),nextstep1:nextstep2)=dnew{i}*gecici*vnew{i}*diag(varR{i})./2-gecici*vnew{i}*diag(driftsR{i});
    gecici=zeros(rownos(i+2)-rownos(i+1));
    for j=1:rownos(i+2)-rownos(i+1)
        gecici(j,j)=exp(dnew{i+1}(j,j)*thresholds(i+1));
    end
    totalmatrix(rownos(i+1)+1:rownos(i+2),nextstep1:nextstep2)=-dnew{i+1}*gecici*vnew{i+1}*diag(varR{i+1})./2+gecici*vnew{i+1}*diag(driftsR{i+1});
end

%I-2
for i=1:r-1
    if it1no(i)>0
        nextstep1=nextstep2+1;
        nextstep2=nextstep2+it1no(i);
        gecici=zeros(rownos(i+1)-rownos(i));
        for j=1:rownos(i+1)-rownos(i)
            gecici(j,j)=exp(dnew{i}(j,j)*thresholds(i+1));
        end
        gecici2=gecici*vnew{i}*diag(varR{i})./2;
        gecici3=zeros(rownos(i+1)-rownos(i),it1no(i));
        for j=1:it1no(i)
            gecici3(:,j)=gecici2(:,it1{i}{j});
        end    
        totalmatrix(rownos(i)+1:rownos(i+1),nextstep1:nextstep2)=gecici3;
        gecici=zeros(rownos(i+2)-rownos(i+1));
        for j=1:rownos(i+2)-rownos(i+1)
            gecici(j,j)=exp(dnew{i+1}(j,j)*thresholds(i+1));
        end
        gecici2=-gecici*vnew{i+1}*diag(varR{i+1})./2;
        gecici3=zeros(rownos(i+2)-rownos(i+1),it1no(i));
        for j=1:it1no(i)
            gecici3(:,j)=gecici2(:,it1{i}{j});
        end            
        totalmatrix(rownos(i+1)+1:rownos(i+2),nextstep1:nextstep2)=gecici3;
    end
end

%I-3
%for i=1:r-1
%    if it2no(i)>0
%        nextstep1=nextstep2+1;
%        nextstep2=nextstep2+it2no(i);
%        gecici=zeros(2*countsplus(i)+countsstar(i));
%        for j=1:2*countsplus(i)+countsstar(i)
%            gecici(j,j)=exp(dnew{i}(j,j)*thresholds(i+1));
%        end
%        gecici2=gecici*vnew{i}*diag(varR{i})/2;
%        gecici3=zeros(rownos(i+1)-rownos(i),it2no(i));
%        for j=1:it2no(i)
%            gecici3(:,j)=gecici2(:,it2{i}{j});
%        end                    
%        totalmatrix(rownos(i)+1:rownos(i+1),nextstep1:nextstep2)=gecici3;
%    end
%end
%I-4
for i=1:r-1
    if it4no(i)+it11no(i)+it13no(i)+it31no(i)+it32no(i)>0
        nextstep1=nextstep2+1;
        nextstep2=nextstep2+it4no(i)+it11no(i)+it13no(i)+it31no(i)+it32no(i);                             
        gecici=zeros(rownos(i+1)-rownos(i));
        for j=1:rownos(i+1)-rownos(i)
            gecici(j,j)=exp(dnew{i}(j,j)*thresholds(i+1));
        end
        gecici2=gecici*vnew{i};
        gecici3=zeros(rownos(i+1)-rownos(i),it4no(i)+it11no(i)+it13no(i)+it31no(i)+it32no(i));
        gecicicount=1;
        for j=1:it11no(i)
            gecici3(:,j)=gecici2(:,it11{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                    
        gecicicount=1;
        for j=it11no(i)+1:it11no(i)+it13no(i)
            gecici3(:,j)=gecici2(:,it13{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                             
        gecicicount=1;
        for j=it11no(i)+it13no(i)+1:it11no(i)+it13no(i)+it31no(i)
            gecici3(:,j)=gecici2(:,it31{i}{gecicicount});
            gecicicount=gecicicount+1;
        end           
        gecicicount=1;
        for j=it11no(i)+it13no(i)+it31no(i)+1:it11no(i)+it13no(i)+it31no(i)+it32no(i)
            gecici3(:,j)=gecici2(:,it32{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                            
        gecicicount=1;
        for j=it11no(i)+it13no(i)+it31no(i)+it32no(i)+1:it11no(i)+it13no(i)+it31no(i)+it32no(i)+it4no(i)
            gecici3(:,j)=gecici2(:,it4{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                            
      
      
        totalmatrix(rownos(i)+1:rownos(i+1),nextstep1:nextstep2)=gecici3;
    end
end

%I-5
%for i=1:r-1
%    if it5no(i)>0
%        nextstep1=nextstep2+1;
%        nextstep2=nextstep2+it5no(i);
%        gecici=zeros(2*countsplus(i+1)+countsstar(i+1));
%        for j=1:2*countsplus(i+1)+countsstar(i+1)
%            gecici(j,j)=exp(dnew{i+1}(j,j)*thresholds(i+1));
%        end
%        gecici2=-gecici*vnew{i+1}*diag(varR{i+1})/2;
%        gecici3=zeros(rownos(i+2)-rownos(i+1),it5no(i));
%        for j=1:it5no(i)
%            gecici3(:,j)=gecici2(:,it2{i}{j});
%        end                    
%        totalmatrix(rownos(i+1)+1:rownos(i+2),nextstep1:nextstep2)=gecici3;
%    end
%end
%I-6

for i=1:r-1
    if it7no(i)+it12no(i)+it13no(i)+it61no(i)+it62no(i)>0
        nextstep1=nextstep2+1;
        nextstep2=nextstep2+it7no(i)+it12no(i)+it13no(i)+it61no(i)+it62no(i);                             
        gecici=zeros(rownos(i+2)-rownos(i+1));
        for j=1:rownos(i+2)-rownos(i+1)
            gecici(j,j)=exp(dnew{i+1}(j,j)*thresholds(i+1));
        end
        gecici2=gecici*vnew{i+1};
        gecici3=zeros(rownos(i+2)-rownos(i+1),it7no(i)+it12no(i)+it13no(i)+it61no(i)+it62no(i));
        gecicicount=1;
        
        for j=1:it12no(i)
            gecici3(:,j)=gecici2(:,it12{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                    
        gecicicount=1;
        for j=it12no(i)+1:it12no(i)+it13no(i)
            gecici3(:,j)=gecici2(:,it13{i}{gecicicount});
            gecicicount=gecicicount+1;
        end           
        gecicicount=1;
        for j=it12no(i)+it13no(i)+1:it12no(i)+it13no(i)+it61no(i)
            gecici3(:,j)=gecici2(:,it61{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                            
        gecicicount=1;
        for j=it12no(i)+it13no(i)+it61no(i)+1:it12no(i)+it13no(i)+it61no(i)+it62no(i)
            gecici3(:,j)=gecici2(:,it62{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                            
        gecicicount=1;
        for j=it12no(i)+it13no(i)+it61no(i)+it62no(i)+1:it12no(i)+it13no(i)+it61no(i)+it62no(i)+it7no(i)
            gecici3(:,j)=gecici2(:,it7{i}{gecicicount});
            gecicicount=gecicicount+1;
        end                                    
        totalmatrix(rownos(i+1)+1:rownos(i+2),nextstep1:nextstep2)=gecici3;
    end
end


%add terms containing C
%T1-I1-T3

%T1
nextstepc1=1;%starting column no
nextstepc2=s;%ending column no
nextstepr1=rownos(r+1);
nextstepr2=rownos(r+1);
if Cat0no>0
    nextstepr1=rownos(r+1)+1;%starting row no
    nextstepr2=rownos(r+1)+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 it1no(i)>0
        nextstepc1=nextstepc2+1;
        nextstepc2=nextstepc2+it1no(i);
    end
end
%I4
for i=1:r-1
    if it4no(i)+it11no(i)+it13no(i)+it31no(i)+it32no(i)>0
        nextstepc1=nextstepc2+1;
        nextstepc2=nextstepc2+it4no(i)+it11no(i)+it13no(i)+it31no(i)+it32no(i);                             
    end
end
%I6
for i=1:r-1
    if it7no(i)+it12no(i)+it13no(i)+it61no(i)+it62no(i)>0
        nextstepc1=nextstepc2+1;
        nextstepc2=nextstepc2+it7no(i)+it12no(i)+it13no(i)+it61no(i)+it62no(i);                             
    end
end

%Normalization conditions
%for C unknowns
totalmatrix(rownos(r+1)+1:rownos(r+1)+Cat0no+sum(Catino),columns)=ones(Cat0no+sum(Catino),1);

%for regimes

for i=1:r-1
    gecici=zeros(rownos(i+1)-rownos(i));
    for j=1:rownos(i+1)-rownos(i)
        if abs(dnew{i}(j,j))<10^-10
            
            gecici(j,j)=thresholds(i+1)-thresholds(i);
        else
            gecici(j,j)=( exp(dnew{i}(j,j)*thresholds(i+1))-exp(dnew{i}(j,j)*thresholds(i)) )/dnew{i}(j,j);
        end
    end
    gecici2=gecici*vnew{i}*ones(s,1);
    totalmatrix(rownos(i)+1:rownos(i+1),columns)=gecici2;
end
gecici=zeros(rownos(r+1)-rownos(r));
for j=1:rownos(r+1)-rownos(r)
    gecici(j,j)=-exp(dnew{r}(j,j)*thresholds(r))/dnew{r}(j,j);        
end
gecici2=gecici*vnew{r}*ones(s,1);
totalmatrix(rownos(r)+1:rownos(r+1),columns)=gecici2;

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


%B)Solution:                                       %|
%-----------                                       %|
totalvector=[zeros(1,columns-1) 1];                    %|
solution=totalvector/totalmatrix;                  %|
%---------------------------------------------------|


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

%pdfs
i=1;                                               %|
for j=1:samples                                    %|
    if j*deltax<=thresholds(i+1)%in regime i                   %|
        gecici=zeros(rownos(i+1)-rownos(i));
        for k=1:rownos(i+1)-rownos(i)
            gecici(k,k)=exp(dnew{i}(k,k)*j*deltax);
        end
        pdfs(:,j)=(solution(rownos(i)+1:rownos(i+1))*gecici*vnew{i})';        
    elseif i<r                                     %|
        i=i+1;
        gecici=zeros(rownos(i+1)-rownos(i));
        for k=1:rownos(i+1)-rownos(i)
            gecici(k,k)=exp(dnew{i}(k,k)*j*deltax);
        end
        pdfs(:,j)=(solution(rownos(i)+1:rownos(i+1))*gecici*vnew{i})';
    end
end                                                %|


%cdfs
i=1;                                               %|
for j=1:samples                                    %|
    if j*deltax<=thresholds(i+1)%in regime i                   %|
        gecici=zeros(rownos(i+1)-rownos(i));
        for k=1:rownos(i+1)-rownos(i)
            if abs(dnew{i}(k,k))<10^-10            
                gecici(k,k)=j*deltax-thresholds(i);
            else
                gecici(k,k)=( exp(dnew{i}(k,k)*j*deltax)-exp(dnew{i}(k,k)*thresholds(i)) )/dnew{i}(k,k);        
            end
        end
        cdfs(:,j)=cdfs(:,j)+(solution(rownos(i)+1:rownos(i+1))*gecici*vnew{i})';        
    elseif i<r                                     %|
        i=i+1;
        for l=j:samples
            cdfs(:,l)=cdfs(:,j-1);
        end
        gecici=zeros(rownos(i+1)-rownos(i));
        for k=1:rownos(i+1)-rownos(i)
            if abs(dnew{i}(k,k))<10^-10            
                gecici(k,k)=j*deltax-thresholds(i);
            else
                gecici(k,k)=( exp(dnew{i}(k,k)*j*deltax)-exp(dnew{i}(k,k)*thresholds(i)) )/dnew{i}(k,k);        
            end
        end
        cdfs(:,j)=cdfs(:,j)+(solution(rownos(i)+1:rownos(i+1))*gecici*vnew{i})';        
    end
end                                                %|

%add accumulations

if Cat0no+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
        
end





S=['b' 'g' 'r' 'c' 'm'];
t=0:deltax:deltax*(samples-1);                     %|
figure(1);                                         %|
hold on                                            %|

for i=1:s                                          %|
    plot(t,pdfs(i,:),S(i));                             %|
end                                                %|
t=0:deltax:deltax*(samples-1);                     %|
figure(2);                                         %|
hold on                                            %|
for i=1:s                                          %|
    plot(t,cdfs(i,:),S(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);                       %|
%---------------------------------------------------|