% Simulation of discrete mousetrap problem
% Input probability of a trap hit, p,
% and number of traps, M
% N = number of balls in the air
% U = unstrung traps
% f = fraction of unsprung traps
p = input(' Probability: p = ');
M = input(' Number of traps: M = ');
f(1) = 1; f(2) = 1 - p/M;
N(1) = 1; U(1) = M;
N(2) = 2*M*(f(1) - f(2));
U(2) = U(1) - N(2)/2;
i = 3;
while ( N(i-1) >= 1 )
    f(i) = f(i-1) + 2*p*f(i-1)*(f(i-1)-f(i-2));
    N(i) = 2*p*N(i-1)*f(i-1);
    U(i) = U(i-1) - N(i)/2;
    i = i + 1;
end
j = i-2; k = i-1;
Bmax = norm(N,inf);
fprintf('Unsprung = %6.2f, bmax = %6.2f, time = %2.0d \n', ...
  U(k),Bmax,j);