function [A0,A1,A2]=examples(examp)

if (examp==1)
    p=10^(-8);
    A0=[1-p 0; 0 0];
    A1=[0 p; 2*p 0];
    A2=[0 0; 0 1-2*p];
end
if (examp==2)
    K=3;

    mks=[2 3 1];
    summks=[2 5 6];

    T{1}=[0 1; 0 0];
    T{2}=[1/2 1/2 0; 0 2/3 1/3; 0 0 1/2];
    T{3}=[4/5];
    Tser=zeros(sum(mks),sum(mks));
    Tser(1:mks(1),1:mks(1))=T{1};
    Tser(mks(1)+1:summks(2),mks(1)+1:summks(2))=T{2};
    Tser(summks(2)+1:sum(mks),summks(2)+1:sum(mks))=T{3};
    
    alpha{1}=[1 0];
    alpha{2}=[1 0 0];
    alpha{3}=[1];

    alphatot=[alpha{1} alpha{2} alpha{3}];
  
    Tser0=ones(sum(mks),1)-sum(Tser')';

    m=3;
    Dund=[999/1000 1/1000 0; 0 999/1000 1/1000; 1/1000 0 999/1000];

    D0=diag([4/5*99/100 99/100*27/28*99/100 99/100*19/20])*Dund;
    D{1}=diag([1/5*99/100 1/100*27/28*99/100 0])*Dund;
    D{2}=diag([0 99/100*1/28*99/100 1/100*19/20])*Dund;
    D{3}=diag([1/100*4/5 99/100*27/28*1/100 99/100*1/20])*Dund;
    D{12}=diag([0 1/100*1/28*99/100 0])*Dund;
    D{13}=diag([1/5*1/100 1/100*27/28*1/100 0])*Dund;
    D{23}=diag([0 99/100*1/28*1/100 1/100*1/20])*Dund;
    D{123}=diag([0 1/100*1/28*1/100 0])*Dund;

    a=4;
    Dtilde0=zeros(a*m,a*m);
    Dtilde0(1:m,1:m)=D0;
    Dtilde{1}=[D{1} D{12} D{13} D{123}; zeros((a-1)*m,a*m)];
    Dtilde{2}=[D{2} zeros(m,m) D{23} zeros(m,m); eye(m) zeros(m,(a-1)*m); ...
           zeros(m,a*m); zeros(m,2*m) eye(m) zeros(m,m)];
    Dtilde{3}=[D{3} zeros(m,(a-1)*m); zeros(m,a*m); eye(m) zeros(m,(a-1)*m); 
        zeros(m,a*m)];

    temp=kron(alpha{1},Dtilde{1}(m+1:end,:));
    for k=2:K
        temp=[temp kron(alpha{k},Dtilde{k}(m+1:end,:))];
    end  
    mtot=size(Tser,1)*m;
    A2=kron(Tser,eye(a*m))+kron(Tser0,[zeros(m,a*mtot); temp]);
    A2=[zeros(m,a*mtot+m); zeros(a*mtot,m) A2];

    temp=kron(alpha{1},Dtilde{1}(1:m,:));
    for k=2:K
        temp=[temp kron(alpha{k},Dtilde{k}(1:m,:))];
    end  
    A1=kron(Tser0,[temp; zeros((a-1)*m,a*mtot)]);
    A1=[zeros(a*mtot+m,m) [temp; A1]];

    A0=[D0 zeros(m,a*mtot); kron(Tser0,[D0; zeros(m*(a-1),m)]) zeros(a*mtot,a*mtot)];
end
if (examp==3)
    k=32;
    delta=10^(-3);
    R=ones(k,k)*(1-delta)/(3*k-3);
    R=R-diag(diag(R));
    A0=R+diag(ones(1,k)*delta);
    A1=R;
    A2=R;
end    
if (examp==4)
    % FS-ALOHA protocol as QBD, arrival rate lambda 
    % T = S+N slots, S for sets, N for resolution
    
    lambda=3;
    D=exp(-lambda);
    som=D;
    i=1;
    while (som < 1 - 10^(-15))
        D(1,i+1)=D(1,i)*lambda/i;
        som=som+D(1,i+1);
        i = i+1;
    end

    load=lambda;

    qm=size(D,2)-2;

    Nslots=5;
    Sslots=5;
    NplusSslots=Nslots+Sslots;


    PN=zeros(qm,qm);
    for q=2:qm+1
        for qac=2:q
            temp=0;
            for v=q-qac:min([q Nslots])
                temp=temp+(-1)^(v+q-qac)*nchoosek(v,q-qac)/Nslots^q*nchoosek(Nslots,v)/factorial(q-v)*...
                    factorial(q)*(Nslots-v)^(q-v);
            end
            PN(q-1,qac-1)=temp;
        end
    end
    PNsolved=zeros(qm,1);
    for q=2:qm+1
        temp=0;
        for v=q:min([q Nslots])
            temp=temp+(-1)^(v+q)*nchoosek(v,q)/Nslots^q*nchoosek(Nslots,v)/factorial(q-v)*...
                factorial(q)*(Nslots-v)^(q-v);
        end
        PNsolved(q-1,1)=temp;
    end

    PS=zeros(qm,qm);
    for q=2:qm+1
        for qac=2:q
            temp=0;
            for v=q-qac:min([q Sslots])
                temp=temp+(-1)^(v+q-qac)*nchoosek(v,q-qac)/Sslots^q*nchoosek(Sslots,v)/factorial(q-v)*...
                    factorial(q)*(Sslots-v)^(q-v);
            end
            PS(q-1,qac-1)=temp;
        end
    end
    PSsolved=zeros(qm,1);
    for q=2:qm+1
        temp=0;
        for v=q:min([q Sslots])
            temp=temp+(-1)^(v+q)*nchoosek(v,q)/Sslots^q*nchoosek(Sslots,v)/factorial(q-v)*...
                factorial(q)*(Sslots-v)^(q-v);
        end
        PSsolved(q-1,1)=temp;
    end

    PNplusS=zeros(qm,qm);
    for q=2:qm+1
        for qac=2:q
            temp=0;
            for v=q-qac:min([q NplusSslots])
                temp=temp+(-1)^(v+q-qac)*nchoosek(v,q-qac)/NplusSslots^q*nchoosek(NplusSslots,v)/factorial(q-v)*...
                    factorial(q)*(NplusSslots-v)^(q-v);
            end
            PNplusS(q-1,qac-1)=temp;
        end
    end
    PNplusSsolved=zeros(qm,1);
    for q=2:qm+1
        temp=0;
        for v=q:min([q NplusSslots])
            temp=temp+(-1)^(v+q)*nchoosek(v,q)/NplusSslots^q*nchoosek(NplusSslots,v)/factorial(q-v)*...
                factorial(q)*(NplusSslots-v)^(q-v);
        end
        PNplusSsolved(q-1,1)=temp;
    end

    pTSgen=1-D(1)-D(2)-D(3:end)*PSsolved;
    TSgen_size=D(3:end)*PS;
    pTSgen_empty=1-D(1)-D(2)-D(3:end)*PNplusSsolved;
    TSgen_size_empty=D(3:end)*PNplusS;
    TSgen_size=TSgen_size/sum(TSgen_size);
    TSgen_size_empty=TSgen_size_empty/sum(TSgen_size_empty);

    % approximate system state at time t

    alpha_ini=1;

    A2=pTSgen*PN;
    A1=pTSgen*PNsolved*TSgen_size+(1-pTSgen)*PN;
    A0=(1-pTSgen)*PNsolved*TSgen_size;
end
if (examp==5)
    %input parameters

    %D-MAP description
    par_c = 10^-3;
    
    lambda1 = 0.6;
    lambda2 = lambda1;
    par_a = 1/50;
    par_b = 1/950-par_c;
    par_d = 6/100;
    par_f = 19*par_c/59;
    par_e = 6/5900-par_f;
    
    D = [   1-par_a         par_a       0           0;
           par_b       1-par_b-par_c   par_c       0;
           0               0           1-par_d     par_d;
           par_f           0           par_e       1-par_e-par_f
       ];
    D1 = diag([lambda1 0 lambda2 0]) * D;
    D0 = D-D1;

    %D = [1-1/500 1/500; 1/1000 1-1/1000];
    %D0 = diag([9/10 3/4])*D;
    %D1 = diag([1/10 1/4])*D;

    %p_input=0.9*0.3;
    %D0=1-p_input;
    %D1=p_input;

    l=size(D,1);

    %service description
    %T=1-0.3;
    %T = [4/5 0 0; 0 1/2 0; 0 0 1/4];
    T = diag(ones(1,9),-1);
    m=size(T,1);
    ml=m*l;
    %alpha=1;
    %alpha = [3/4 1/8 1/8];
    alpha = [2.6012295*10^(-2)  1.1705533*10^(-1)  2.3411065*10^(-1)  2.7312910*10^(-1)  2.0484682*10^(-1)  1.0242341*10^(-1)  3.4141137*10^(-2) 7.3159579*10^(-3)  9.1449474*10^(-4)  5.0805263*10^(-5)];
    alpha=alpha/sum(alpha); % this input alfa was not 100% exact, hence the normalization
    T0 = ones(m,1)-T*ones(m,1);

    % approximate queuelength the time t

    A2=kron(D1,T);
    A1=kron(D0,T)+kron(D1,T0*alpha);
    A0=kron(D0,T0*alpha);
end