format compact;

for mu = 2:.5:4.5;
mu
r = zeros(5,5);   
p = zeros(9,1);

fct = 1;

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

    if i >=1;    
    p(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 p]

r(1,1) = p(1,1)+p(2,1)+p(3,1)+p(4,1);     % p00 %
r(1,2) = p(5,1);                          % p01 %
r(1,3) = p(6,1);                          % p02 %
r(1,4) = p(7,1);                          % p03 %
r(1,5) = p(8,1);                          % p04 %

r(2,1) = p(1,1)+p(2,1)+p(3,1);            % p10 %
r(2,2) = p(4,1);                          % p11 %
r(2,3) = p(5,1);                          % p12 %
r(2,4) = p(6,1);                          % p13 %
r(2,5) = p(7,1);                          % p14 %

r(3,1) = p(1,1)+p(2,1);                   % p20 %
r(3,2) = p(3,1);                          % p21 %
r(3,3) = p(4,1);                          % p22 %
r(3,4) = p(5,1);                          % p23 %
r(3,5) = p(6,1);                          % p24 %

r(4,1) = p(1,1);                          % p30 %
r(4,2) = p(2,1);                          % p31 %
r(4,3) = p(3,1);                          % p32 %
r(4,4) = p(4,1);                          % p33 %
r(4,5) = p(5,1);                          % p34 %

r(5,1) = 0;                               % p40 %
r(5,2) = p(1,1);                          % p41 %
r(5,3) = p(2,1);                          % p42 %
r(5,4) = p(3,1);                          % p43 %
r(5,5) = p(4,1);                          % p44 %
r

ww= eye(5) - r;
w = inv(ww)

arlvec = sum(w');
x0 = [0;1;2;3;4];
arl = [x0 arlvec']
end;
end