function yp = mouseODE(t,y,p,M)
% Mousetrap ODE (with M traps)
C = (2/p)*log(2*p+1)-2-p/M;
yp = 2*y - (2/p)*log(2*p*y + 1) + C;
end