clear
K = 75;

pe = 0.99;
B = [ 1 zeros(1,K-2)];
S0 = [ zeros(K-2,1); pe ]
S = zeros(K-1,K-1);

for i=1:K-1
    S(i,i) = 1-pe;
    if (i < K-1)        
        S(i,i+1) = pe;
    end
end
S
total = 0;
i = K;
go = 1;
while ((total < 0.999) && (i < 500) && (go == 1))
    calc = choose(i-1,K-1)*(1-pe)^(i-K)*pe^(K);
            
    if ((calc ~= NaN) && (calc ~= Inf))
        total = total + calc;
        go = 1;
        Dist(i-K+1) = calc;
        Num(i-K+1) = i;
        if ((i-K+1) == 1)
            Cumlative(i-K+1) = 0;
        else
            Cumlative(i-K+1) = Cumlative(i-K) + calc;
        end
    else
        go = 0;        
        i        
    end;
    i = i+1;
end

bar(Num,Dist)
%figure
%bar(Num,Cumlative)

S_prime = [ 0 B;
            S0 S ]

[ order_S_prime, null ] = size(S_prime);

E = zeros(order_S_prime*order_S_prime,order_S_prime*order_S_prime+1); %initialze the matrix holding the set of equations 
for i=1:order_S_prime
    for j=1:order_S_prime
        E((i-1)*order_S_prime+j,(i-1)*order_S_prime+j) = 1;
        for k=1:order_S_prime
            if (k ~= j)
                % this is the mean first recurrence calculation
                E((i-1)*order_S_prime+j,(k-1)*order_S_prime+j) = E((i-1)*order_S_prime+j,(k-1)*order_S_prime+j) - S_prime(i,k);
            end
        end
    end     
end
for h=1:order_S_prime*order_S_prime    
    E(h,(order_S_prime*order_S_prime)+1) = 1; % add ones to the last column
end

E= rref(E); %solve the equations using Gaussian elimination

for i=1:order_S_prime
    for j=1:order_S_prime
        M(i,j) = E((i-1)*order_S_prime+j,order_S_prime*order_S_prime+1); %arrange the results in a matrix
    end
end

mean_time_to_absorbtion = M(1,1)

pi = [ 1 zeros(1,order_S_prime-1) ];

error = zeros(1,order_S_prime);
loop = 1; %this will be set to 1 if the error of any term is greater than 10^-12
while (loop == 1)
    pi_prime = pi*S_prime;
    error = pi_prime-pi;
    pi = pi_prime;  
    loop = 0;
    for i=1:order_S_prime
        if (abs(error(i)) > 10^-12)
            loop = 1;
        end               
    end
end

pi
