clear
%---------------------------------------------------|
%----MULTI REGIME - FIRST ORDER - FINITE BUFFER-----|
%---------WITH ORDERED SCHUR DECOMPOSITION----------|
%-------------THE MOST GENERAL CASE-----------------|
%-------------------INCLUDING-----------------------|
%-------------REPULSIVE TYPE BOUNDARIES-------------|
%--------------&ZERO RATE EXTENSIONS----------------|
%---------------------------------------------------|

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

clear
%------------------PARAMETERS-----------------------|                                                   

%***************************EXAMPLE 1************************
Thresholds=[0,1,2]; %terminal and intermediate boundary points
ratesR{1}=[1,-2,2,-4,0];%rate in region 1: states 1 and 2 respectively     
ratesR{2}=[-1,1,3,-1,2];%rate in 
RT{1}=[1,0,2,0,0];%rate in region 1: states 1 and 2 respectively     
RT{2}=[0,0,2,-4,0];%rate in 
RT{3}=[-1,0,0,-1,0];%rate in region 1: states 1 and 2 respectively  

QR{1}=[-4,1,1,1,1;1,-4,1,1,1;1,1,-4,1,1;1,1,1,-4,1;1,1,1,1,-4];%infinitesimal generator in region 1          
QR{2}=QR{1};%[-5,2,1,1,1;0,-3,1,1,1;2,1,-5,0,2;0,1,1,-3,1;1,1,1,1,-4];%infinitesimal generator in region 1                    
QT{1}=QR{1};%infinitesimal generator in region 1          
QT{2}=QR{1};%infinitesimal generator in region 2          
QT{3}=QR{2};%infinitesimal generator in region 1          
%***************************EXAMPLE 2************************

boundarytype=1;                                    %| 
%---------------------------------------------------|


%-------------------PREPARE-------------------------|
r=length(ratesR)                                  %|
s=length(ratesR{1})                               %|
%A)Find States of Zeros&Number of Zeros:           %|
%---------------------------------------           %|
%zeroR{i}(j): state of jth zero at region i        %|
for i=1:r                                          %|
                                           %|
    zeroR{i}(1)=-1;                         %|
    k=2;
    for j=1:s                                      %|
        if ratesR{i}(j)==0                         %|
            zeroR{i}(k)=j;                         %|
            k=k+1;                                 %|
        end                                        %|
    end                                            %|
    
end  
zeroR%|
%noofzeros(i): number of zeros in ith region       %|
noofzeros=zeros(1,r);                              %|
for i=1:r  
    %|
    noofzeros(i)=length(zeroR{i})-1;                 %|
end                                                %|
%noofzeros2(i):number of zeros upto boundary i     %|

noofzeros2=zeros(1,r+1);                           %|
for i=1:r                                          %|
    noofzeros2(i+1)=sum(noofzeros(1:i));           %|
end
noofzeros
s
r
noofzeros2

%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%B)Find Rates Without Zeros:                       %|
%---------------------------                       %|
%ratesnz{i}=rates in region i without zeros        %|              
for i=1:r                                          %|
    ratesnz{i}=zeros(1,s-noofzeros(i));            %|
    k=1;                                           %|
    for j=1:s                                      %|
        if ~(ratesR{i}(j)==0)                      %|
            ratesnz{i}(k)=ratesR{i}(j);            %|
            k=k+1;                                 %|
        end                                        %|
    end                                            %|
end                                                %|
ratesnz{1}
ratesnz{2}
%ratesnz{3}
%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%C)Find Rates&Q matrices Without Zeros:            %|
%--------------------------------------            %|
%Q=[Qnz Qa;Q0 Qb]                                  %|
%Qnew{i}:infinitesimal generator of nonzero states %|
for rr=1:r                                         %|
    R=diag(ratesnz{rr});                           %| 
    Qb=zeros(noofzeros(rr));                       %|
    Qa=zeros(s-noofzeros(rr),noofzeros(rr));       %|
    Q0=zeros(noofzeros(rr),s-noofzeros(rr));       %|
    Qnz=zeros(s-noofzeros(rr),s-noofzeros(rr));    %|
    ib=1;                                          %|
    ia=1;                                          %|
    i0=1;                                          %|
    inz=1;                                         %|
    for i=1:s                                      %|
        jb=1;                                      %|
        ja=1;                                      %|
        j0=1;                                      %|
        jnz=1;                                     %|
        for j=1:s                                  %|   
            if ratesR{rr}(i)==0&ratesR{rr}(j)==0   %|
                Qb(ib,jb)=QR{rr}(i,j);             %|
                jb=jb+1;                           %|
            elseif ~(ratesR{rr}(i)==0)&ratesR{rr}(j)==0
                Qa(ia,ja)=QR{rr}(i,j);             %|
                ja=ja+1;                           %|
            elseif ratesR{rr}(i)==0&~(ratesR{rr}(j)==0)
                Q0(i0,j0)=QR{rr}(i,j);             %|
                j0=j0+1;                           %|
            else                                   %|
                Qnz(inz,jnz)=QR{rr}(i,j);          %|
                jnz=jnz+1;                         %|
            end                                    %|
        end                                        %|
        if ratesR{rr}(i)==0                        %|
            ib=ib+1;                               %|
            i0=i0+1;                               %|
        else                                       %|
            ia=ia+1;                               %|
            inz=inz+1;                             %|
        end                                        %|
    end                                            %|
    QnewR{rr}=Qnz-Qa*(Qb)^-1*Q0;                   %|
    %carpim: fnz*carpim=f0                         %|
    carpim{rr}=-Qa*Qb^-1;                          %|
end                                                %|
%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%C)Find potentially nonzero accumulations:         %|
%-----------------------------------------         %|
%noofpnz: number of potentially nonzero boundaries %|
%pnzstate2: vector of states of pnz bound.         %|
%pnzbound2: vector of boundaries of pnz bound.     %|
noofpnz=0;                                         %|
pnzstate=zeros(1,s*r+s);                           %|
pnzbound=zeros(1,s*r+s);                           %|
for i=1:s                                          %|
    for j=1:r                                      %|
        if ratesR{j}(i)==0                         %|
            if j==1                                %|
                if ratesR{2}(i)>0                  %|
                    noofpnz=noofpnz+1;             %|
                    pnzbound(noofpnz)=j;           %|
                    pnzstate(noofpnz)=i;           %|
                end                                %|
            elseif j==r                            %|
                if ratesR{r-1}(i)<0                %|
                    noofpnz=noofpnz+1;             %|
                    pnzbound(noofpnz)=j-1;         %|
                    pnzstate(noofpnz)=i;           %|
                end                                %|
            else                                   %|
                if ratesR{j+1}(i)>0                %|
                    noofpnz=noofpnz+1;             %|
                    pnzbound(noofpnz)=j;           %|
                    pnzstate(noofpnz)=i;           %|
                end                                %|
                if ratesR{j-1}(i)<0                %|
                    noofpnz=noofpnz+1;             %|
                    pnzbound(noofpnz)=j-1;         %|
                    pnzstate(noofpnz)=i;           %|
                end                                %| 
            end                                    %|
        end                                        %|
    end                                            %|
end                                                %|
pnzstate2=pnzstate(1:noofpnz);                     %|

pnzbound2=pnzbound(1:noofpnz);                     %|
pnzstate2
pnzbound2
%D)Find Repulsive Boundaries:                      %|
%----------------------------                      %|
rep=0;                                             %| 
if r>1                                             %|
    for i=1:s                                      %|
        for j=1:r-1                                %|
            if ratesR{j}(i)<0&ratesR{j+1}(i)>0     %|
                rep=rep+1;                         %|
                repstate{rep}=i;                   %|                
                repregion{rep}=j;                  %|  
            end                                    %| 
        end                                        %|
    end                                            %|
end                                                %|
rep

%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%E)Change Number of Variables Acc. to Boundaries:  %|
%------------------------------------------------  %|
addvar=rep;                                        %|
if boundarytype==1                                 %|
    addvar=2*rep;                                  %|
end                                                %|
accsbound=zeros(1,s+addvar+noofpnz+sum(noofzeros));%takesvalues1->Nregions
accsstate=zeros(1,s+addvar+noofpnz+sum(noofzeros));%takesvalues1->Nstates 
%accno: number of accumulations                    %|
accno=s+addvar+noofpnz+sum(noofzeros);             %|
%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%F)Find Accumulations:                             %|
%---------------------                             %|
k=1;                                               %|
%Accumulations at zero                             %|
for i=1:s                                          %|
    if ratesR{1}(i)<0                              %|
        accsstate(k)=i;                            %|
        accsbound(k)=0;                            %|
        k=k+1;                                     %|
    end                                            %|
end      
accsbound%|
accsstate
%Accumulation at Intermediate Boundaries           %|
for j=1:r-1                                        %|
    for i=1:s                                      %|
        if ratesR{j}(i)>0&ratesR{j+1}(i)<0         %|
            accsstate(k)=i;                        %|
            accsbound(k)=j;                        %|
            k=k+1;                                 %|
        end                                        %|
    end                                            %|
end                                                %|
if boundarytype==1                                 %|
    for j=1:r-1                                    %|
        for i=1:s                                  %|
            if ratesR{j}(i)<0&ratesR{j+1}(i)>0     %|
                accsstate(k)=i;                    %|
                accsbound(k)=j;                    %|
                k=k+1;                             %|
            end                                    %|
        end                                        %|
    end                                            %|
end                                                %|
accsbound%|
accsstate

%Accumulations at Buffer Length                    %|
for i=1:s                                          %|
    if ratesR{r}(i)>0                              %|
        accsstate(k)=i;                            %|
        accsbound(k)=r;                            %|
        k=k+1;                                     %|
    end                                            %|
end
ratesR{1}

ratesR{2}
%ratesR{3}
accsbound%|
accsstate
%|
%Accumulations at Both sides of regions with zero rate
for i=1:r                                          %|
    for j=1:s                                      %|
        if ratesR{i}(j)==0                         %|
            accsstate(k)=j;                        %|
            accsbound(k)=i;                        %|
            k=k+1;                                 %|
            if i>1                                 %|
                if ~(ratesR{i-1}(j)==0)            %|
                    accsstate(k)=j;                %|
                    accsbound(k)=i-1;              %|
                    k=k+1;                         %|
                end                                %|
            else                                   %|
                accsstate(k)=j;                    %|
                accsbound(k)=i-1;                  %|
                k=k+1;                             %|
            end                                    %|
        end                                        %|
    end                                            %|
end                                                %|
accsbound%|
accsstate

%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%G)Ordered Schur Decomposition:                    %|
%------------------------------                    %|
for i=1:r                                          %|
   %[U,U_inv,A,Aplus,Aminus,Nminus,pik,lambda]=additive2_decompose_schur2(diag(ratesnz{i}),QnewR{i});
   [U_inv,Aplus,Aminus,Nminus] = additive2_decompose_schur(diag(ratesnz{i}),QnewR{i});
   Nplus=(s-noofzeros(i))-Nminus-1;                %|
   Lzero=U_inv(Nminus+1,1:(s-noofzeros(i)));       %|
   Lminus=U_inv(1:Nminus,1:(s-noofzeros(i)));      %|
   Lplus=U_inv(Nminus+2:(s-noofzeros(i)),1:(s-noofzeros(i)));                    %|
   Uinvc{i}=U_inv;                                 %|
   Aplusc{i}=Aplus;                                %|
   Aminusc{i}=Aminus;                              %|
   Nplusc{i}=Nplus;                                %|
   Nminusc{i}=Nminus;                              %|
   Lzeroc{i}=Lzero;                                %|
   Lplusc{i}=Lplus;                                %|
   Lminusc{i}=Lminus;                              %|
end                                                %|
%---------------------------------------------------|

sutun=sum(noofzeros);
%---------------------SOLUTION----------------------|
totalmatrix=zeros(r*s-sutun+accno,r*s+s+1+addvar+noofpnz);   
size(totalmatrix)
%A) Boundary Conditions:                           %|
%-----------------------                           %|
%Boundary C. at Zero                               %| 
a=[Lzeroc{1}*diag(ratesnz{1});...                  %|
 Lminusc{1}*diag(ratesnz{1});...                   %|
 expm(-Aplusc{1}*Thresholds(2))*Lplusc{1}*diag(ratesnz{1})]; 
a2=zeros(s-noofzeros(1),s);                        %|
if noofzeros(1)>0                                  %|
    k=1;                                           %|
    m=2;                                           %|
    for ii=1:s                                     %|
        
        if ~(zeroR{1}(m)==ii)                      %|
            a2(:,ii)=a(:,k);                       %|
            k=k+1;                                 %|
        else                                       %|
             if noofzeros(1)>m                     %|
                m=m+1;                             %|
             end                                   %|
        end                                        %|
    end                                            %|
else                                               %|
    a2=a;                                          %|
end                                                %|
totalmatrix(1:(s-noofzeros(1)),1:s)=a2;            %|
integ1=[Lzeroc{1}*Thresholds(2);...                %|
 Aminusc{1}^-1*(expm(Aminusc{1}*Thresholds(2))-eye(Nminusc{1}))*Lminusc{1};...
 Aplusc{1}^-1*(eye(Nplusc{1})-expm(-Aplusc{1}*Thresholds(2)))*Lplusc{1}];
if noofzeros(1)>0
integ2=integ1*ones(s-noofzeros(1),1)+integ1*carpim{1}*ones(noofzeros(1),1);
else
    integ2=integ1*ones(s-noofzeros(1),1);
end
totalmatrix(1:s-noofzeros(1),r*s+s+1)=integ2;      %|
totalmatrix
%Boundary C. at Thresholds                         %|
if r>1                                             %|
    for i=1:r-1                                    %|
        a=[-Lzeroc{i}*diag(ratesnz{i});...         %|
         -expm(Aminusc{i}*(Thresholds(i+1)-Thresholds(i)))*Lminusc{i}*diag(ratesnz{i});...
         -Lplusc{i}*diag(ratesnz{i})];             %|
        a2=zeros(s-noofzeros(i),s);                %|
        if noofzeros(i)>0                          %|
            k=1;                                   %|
            m=2;                                   %|
            for ii=1:s                             %|
                if ~(zeroR{i}(m)==ii)              %|
                    a2(:,ii)=a(:,k);               %|
                    k=k+1;                         %|
                else                               %|
                    if noofzeros(i)>m              %|
                        m=m+1;                     %|
                    end                            %|
                end                                %|
            end                                    %|
        else                                       %|
            a2=a;                                  %|
        end                                        %|
        totalmatrix((i-1)*s-noofzeros2(i)+1:i*s-noofzeros2(i+1),i*s+1:(i+1)*s)=a2;
        a=[Lzeroc{i+1}*diag(ratesnz{i+1});...      %|
         Lminusc{i+1}*diag(ratesnz{i+1});...       %|
         expm(-Aplusc{i+1}*(Thresholds(i+2)-Thresholds(i+1)))*Lplusc{i+1}*diag(ratesnz{i+1})];
        a2=zeros(s-noofzeros(i+1),s);              %|
        if noofzeros(i+1)>0                        %|
            k=1;                                   %|
            m=2;                                   %|
            for ii=1:s                             %|
                if ~(zeroR{i+1}(m)==ii)            %|
                    a2(:,ii)=a(:,k);               %|
                    k=k+1;                         %|
                else                               %|
                    if noofzeros(i+1)>m            %|
                        m=m+1;                     %| 
                    end                            %|
                end                                %|
            end                                    %|
        else                                       %|
            a2=a;                                  %|
        end                                        %|
        totalmatrix(i*s+1-noofzeros2(i+1):(i+1)*s-noofzeros2(i+2),i*s+1:(i+1)*s)=a2;
        integ1=[Lzeroc{i+1}*(Thresholds(i+2)-Thresholds(i+1));...
         Aminusc{i+1}^-1*(expm(Aminusc{i+1}*(Thresholds(i+2)-Thresholds(i+1)))-eye(Nminusc{i+1}))*Lminusc{i+1};...
         Aplusc{i+1}^-1*(eye(Nplusc{i+1})-expm(-Aplusc{i+1}*(Thresholds(i+2)-Thresholds(i+1))))*Lplusc{i+1}];
 disp('1')       
 integ1*ones(s-noofzeros(i+1),1)
 disp('2')
 integ1*carpim{i+1}*ones(noofzeros(i+1))
 if noofzeros(i+1)>0
     
        integ2=integ1*ones(s-noofzeros(i+1),1)+integ1*carpim{i+1}*ones(noofzeros(i+1));
    else
        integ2=integ1*ones(s-noofzeros(i+1),1);
    end
        totalmatrix(i*s+1-noofzeros2(i+1):(i+1)*s-noofzeros2(i+2),r*s+s+1)=integ2;  %|
    end                                            %|
else                                               %|
    i=0;                                           %|
end        
totalmatrix%|
%Boundary C. at Buffer Length                      %| 
a=[Lzeroc{i+1}*diag(ratesnz{i+1});...              %|
 expm(Aminusc{i+1}*(Thresholds(i+2)-Thresholds(i+1)))*Lminusc{i+1}*diag(ratesnz{i+1});...
 Lplusc{i+1}*diag(ratesnz{i+1})];                  %|
a2=zeros(s-noofzeros(i+1),s);                      %|
if noofzeros(i+1)>0                                %|
    k=1;                                           %|
    m=2;                                           %| 
    for ii=1:s                                     %|
       if ~(zeroR{i+1}(m)==ii)                     %|
            a2(:,ii)=a(:,k);                       %|
            k=k+1;                                 %|
       else                                        %|
            if noofzeros(i+1)>m                    %|
                m=m+1;                             %|
            end                                    %|
       end                                         %|
    end                                            %| 
else                                               %|
    a2=a;                                          %|
end                                                %|   
totalmatrix(i*s+1-noofzeros2(i+1):(i+1)*s-noofzeros2(i+2),(i+1)*s+1:(i+2)*s)=a2;
totalmatrix((i+1)*s-noofzeros2(i+2)+1:(i+1)*s-noofzeros2(i+2)+accno,r*s+s+1)=ones(accno,1);      %|
pnzbound2
pnzstate2
accsbound
accsstate
if noofpnz>0                                       %|                                        
    for j=1:noofpnz                                %|
        for k=1:accno                              %| 
            if accsbound(k)==pnzbound2(j)          %|
                accsstate(k)
                pnzstate2(j)
                totalmatrix((i+1)*s-noofzeros2(i+2)+k,r*s+s+1+j+addvar)=QT{pnzbound2(j)+1}(accsstate(k),pnzstate2(j));
                QT{pnzbound2(j)+1}
            end                                    %|
        end                                        %|  
    end                                            %|
end

%totalmatrix%|
for i=1:accno                                      %|
    if accsbound(i)==0                             %|
        totalmatrix(r*s-sutun+i,1:s)=-QT{1}(accsstate(i),:);
    elseif accsbound(i)==r                         %|
        totalmatrix(r*s-sutun+i,r*s+1:(r+1)*s)=QT{r+1}(accsstate(i),:);     
    else                                           %|
        totalmatrix(r*s-sutun+i,accsbound(i)*s+1:(accsbound(i)+1)*s)=-QT{accsbound(i)+1}(accsstate(i),:);
    end                                            %|
end   
accsbound
accsstate%|
totalmatrix

%totalmatrix
%Boundary C. For Repulsive Boundaries              %|
if rep>0                                           %|
%    for j=1:rep                                    %|
%        i=repregion{j};                            %|
%        X=[Lzeroc{i}*diag(ratesR{i});...           %|
%         expm(Aminusc{i}*(Thresholds(i+1)-Thresholds(i)))*Lminusc{i}*diag(ratesR{i});...
%         Lplusc{i}*diag(ratesR{i})];               %|
%        X2=X(:,repstate{j});                       %|
%        if boundarytype==1                         %|
%            totalmatrix((i-1)*s+1:i*s,r*s+s+(j-1)*2+2)=X2;
%        elseif boundarytype==2                     %|
%            totalmatrix((i-1)*s+1:i*s,r*s+s+(j-1)+2)=X2;
%        end                                        %|
%        X=[Lzeroc{i+1}*diag(ratesR{i+1});...       %|
%         Lminusc{i+1}*diag(ratesR{i+1});...        %|
%         expm(-Aplusc{i+1}*(Thresholds(i+2)-Thresholds(i+1)))*Lplusc{i+1}*diag(ratesR{i+1})];
%        X2=X(:,repstate{j});                       %|
%        if boundarytype==1                         %|
%            totalmatrix(i*s+1:(i+1)*s,r*s+s+2*j+1)=X2;
%        elseif boundarytype==3                     %|
%            totalmatrix(i*s+1:(i+1)*s,r*s+s+j+1)=X2;
%        end                                        %|
%    end                                            %|
    for j=1:rep
        i=repregion{j};
        k=repstate{j};
        totalmatrix((i-1)*s-noofzeros2(i)+1:i*s-noofzeros2(i+1),r*s+s+j*2)=totalmatrix((i-1)*s-noofzeros2(i)+1:i*s-noofzeros2(i+1),i*s+k);
        i=i+1;
        totalmatrix((i-1)*s-noofzeros2(i)+1:i*s-noofzeros2(i+1),r*s+s+j*2+1)=totalmatrix((i-1)*s-noofzeros2(i)+1:i*s-noofzeros2(i+1),(i-1)*s+k);
    end
        
        
end                                                %|
totalmatrix
%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%totalmatrix

%B)Solution:                                       %|
%-----------                                       %|
totalvector=[zeros(1,r*s+s) 1 zeros(1,addvar) zeros(1,noofpnz)];    %|
solution=totalvector/totalmatrix;                  %|
%---------------------------------------------------|


solution
%---------------------PLOTS-------------------------|
samples=9000+1;                                    %|
deltax=Thresholds(r+1)/9000;                       %|
pdfs=zeros(s,samples);                             %|
cdfs=zeros(s,samples);                             %|
totalpdfs=zeros(1,samples);                        %|
totalcdfs=zeros(1,samples);                        %|
for i=1:r                                          %|
    azeroc{i}=solution((i-1)*s-noofzeros2(i)+1);                 %|
    aminusc{i}=solution((i-1)*s-noofzeros2(i)+2:(i-1)*s-noofzeros2(i)+1+Nminusc{i});
    aplusc{i}=solution((i-1)*s-noofzeros2(i)+2+Nminusc{i}:i*s-noofzeros2(i+1));  %|
end                                                %|
i=1;                                               %|
for j=1:samples                                    %|
    if j*deltax<=Thresholds(i+1)                   %|
            a=(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});            %|
            b=a*carpim{i};
            a2=zeros(1,s);                %|
            if noofzeros(i)>0                          %|
                k=1;                                   %|
                m=1;                                   %|
                for ii=1:s                             %|
                    if ~(zeroR{i}(m+1)==ii)              %|
                        a2(ii)=a(k);               %|
                        k=k+1;                         %|
                    else
                        a2(ii)=b(m);               %|%|
                        if noofzeros(i)>m              %|
                            m=m+1;                     %|
                        end                            %|
                    end                                %|
                end                                    %|
            else                                       %|
            a2=a;                                  %|
            end                                        %|

        
        
        pdfs(:,j)=a2';
    elseif i<r                                     %|
        i=i+1;                                     %|
            a=(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});            %|
            b=a*carpim{i};
            a2=zeros(1,s);                %|
            if noofzeros(i)>0                          %|
                k=1;                                   %|
                m=1;                                   %|
                for ii=1:s                             %|
                    if ~(zeroR{i}(m+1)==ii)              %|
                        a2(ii)=a(k);               %|
                        k=k+1;                         %|
                    else
                        a2(ii)=b(m);               %|%|
                        if noofzeros(i)>m              %|
                            m=m+1;                     %|
                        end                            %|
                    end                                %|
                end                                    %|
            else                                       %|
            a2=a;                                  %|
            end                                        %|


        
        
        pdfs(:,j)=a2';
    else                                           %|

            a=(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});            %|
            b=a*carpim{i};
            a2=zeros(1,s);                %|
            if noofzeros(i)>0                          %|
                k=1;                                   %|
                m=1;                                   %|
                for ii=1:s                             %|
                    if ~(zeroR{i}(m+1)==ii)              %|
                        a2(ii)=a(k);               %|
                        k=k+1;                         %|
                    else
                        a2(ii)=b(m);               %|%|
                        if noofzeros(i)>m              %|
                            m=m+1;                     %|
                        end                            %|
                    end                                %|
                end                                    %|
            else                                       %|
            a2=a;                                  %|
            end                                        %|


        
        
        pdfs(:,j)=a2';
    end                                            %|
    
end                                                %|
t=0:deltax:deltax*(samples-1);                     %|
figure(1);                                         %|
hold on                                            %|
for i=1:s                                          %|
    plot(t,pdfs(i,:));                             %|
end                                                %|
totalpdfs=ones(1,s)*pdfs;                          %|
t=0:deltax:deltax*(samples-1);                     %|
figure(2);                                         %|
hold on                                            %|
plot(t,totalpdfs);                                 %|

i=1;                                               %|
mm=0;                                               %|
for j=1:samples                                    %|
    if j*deltax<=Thresholds(i+1)                   %|
                    a=(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});
            b=a*carpim{i};
            a2=zeros(1,s);                %|
            if noofzeros(i)>0                          %|
                k=1;                                   %|
                m=1;                                   %|
                for ii=1:s                             %|
                    if ~(zeroR{i}(m+1)==ii)              %|
                        a2(ii)=a(k);               %|
                        k=k+1;                         %|
                    else
                        a2(ii)=b(m);               %|%|
                        if noofzeros(i)>m              %|
                            m=m+1;                     %|
                        end                            %|
                    end                                %|
                end                                    %|
            else                                       %|
            a2=a;                                  %|
            end                                        %|

        

        cdfs(:,j)= mm+a2';  
    elseif i<r                                     %|
        mm=cdfs(:,j-1);                             %|
        i=i+1;                                     %|
       a=(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});
            b=a*carpim{i};
            a2=zeros(1,s);                %|
            if noofzeros(i)>0                          %|
                k=1;                                   %|
                m=1;                                   %|
                for ii=1:s                             %|
                    if ~(zeroR{i}(m+1)==ii)              %|
                        a2(ii)=a(k);               %|
                        k=k+1;                         %|
                    else
                        a2(ii)=b(m);               %|%|
                        if noofzeros(i)>m              %|
                            m=m+1;                     %|
                        end                            %|
                    end                                %|
                end                                    %|
            else                                       %|
            a2=a;                                  %|
            end                                        %|


        cdfs(:,j)= mm+a2';
    else
                            a=(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});
            b=a*carpim{i};
            a2=zeros(1,s);                %|
            if noofzeros(i)>0                          %|
                k=1;                                   %|
                m=1;                                   %|
                for ii=1:s                             %|
                    if ~(zeroR{i}(m+1)==ii)              %|
                        a2(ii)=a(k);               %|
                        k=k+1;                         %|
                    else
                        a2(ii)=b(m);               %|%|
                        if noofzeros(i)>m              %|
                            m=m+1;                     %|
                        end                            %|
                    end                                %|
                end                                    %|
            else                                       %|
            a2=a;                                  %|
            end                                        %|


        cdfs(:,j)= mm+a2';
    end                                            %|
end                                                %|   
accbef=0;                                          %|
for k=1:accno                                          %|
   if accsbound(k)==0                              %|
        accbef=solution(s*r-noofzeros2(r+1)+k);                    %|
        cdfs(accsstate(k),1:samples)=cdfs(accsstate(k),1:samples)+accbef;
   end                                             %| 
end                                                %|
for i=1:r-1                                        %|
    for k=1:accno                                      %|
        if accsbound(k)==i                         %|
            accbef=solution(s*r-noofzeros2(r+1)+k);                %|
            cdfs(accsstate(k),floor(Thresholds(i+1)/deltax):samples)=cdfs(accsstate(k),floor(Thresholds(i+1)/deltax):samples)+accbef;    
        end                                        %|
    end                                            %|
end                                                %|
for k=1:accno                                          %|
    if accsbound(k)==r                             %|
        accbef=solution(s*r-noofzeros2(r+1)+k);                    %|
        cdfs(accsstate(k),samples)=cdfs(accsstate(k),samples-1)+accbef;
    end                                            %|
end                                                %|
t=0:deltax:deltax*(samples-1);                     %|
figure(3);                                         %|
hold on                                            %|
for i=1:s                                          %|
    plot(t,cdfs(i,:),'y');                             %|
end                                                %|
totalcdfs=ones(1,s)*cdfs;                          %|
t=0:deltax:deltax*(samples-1);                     %|
hold on                                            %|
figure(4);hold on;plot(t,totalcdfs);                       %|
%---------------------------------------------------|
