clear all
close all
n = 5;
J = 10000;
P = catenanm(n);
N = 10;
i = 3;
v = zeros(1,n);
v(i) = 1;
for k = 1:N
  v = v*P;
end
f = zeros(1,n);
for j = 1:J
  k = simula(P,N,i);
  f(k) = f(k)+1;
end
f = f/J;
bar([1:n],f)
hold on
plot([1:n],v,'r*')