format compact;

for mu = 2:1:4;
mu
pr = zeros(9,1);
p  = zeros(6,6);

p(6,6) = 1;
for i = 1:5; p(6,i) = 0; end;
           
fct = 1;

for i = 0:8;
    if i ==0;    
    pr(i+1,1) = exp(-mu)*(mu^i);
    end;

    if i >=1;    
    pr(i+1,1) = exp(-mu)*(mu^i)/fct;
    fct = fct*(i+1);
    end;
end;
x=[0;1;2;3;4;5;6;7;8];
dist = [x pr]

p(1,1) = pr(1,1)+pr(2,1)+pr(3,1)+pr(4,1);
p(1,2) = pr(5,1);
p(1,3) = pr(6,1);
p(1,4) = pr(7,1);
p(1,5) = pr(8,1);

p(2,1) = pr(1,1)+pr(2,1)+pr(3,1);
p(2,2) = pr(4,1);
p(2,3) = pr(5,1);
p(2,4) = pr(6,1);
p(2,5) = pr(7,1);

p(3,1) = pr(1,1)+pr(2,1);
p(3,2) = pr(3,1);
p(3,3) = pr(4,1);
p(3,4) = pr(5,1);
p(3,5) = pr(6,1);

p(4,1) = pr(1,1);
p(4,2) = pr(2,1);
p(4,3) = pr(3,1);
p(4,4) = pr(4,1);
p(4,5) = pr(5,1);

p(5,1) = 0;      
p(5,2) = pr(1,1);
p(5,3) = pr(2,1);
p(5,4) = pr(3,1);
p(5,5) = pr(4,1);

for i = 1:5; 
    rsum = 0;
for j = 1:5; 
    rsum = rsum + p(i,j);
end;
    p(i,6) = 1- rsum;
end;

p
p2 = p*p
p4 = p2*p2
p8 = p4*p4
p16 = p8*p8
p32 = p16*p16
end;
end;
end