It is almost always a bad idea to generate random variates one at a time in a time-sensitive construct in MM. This is a rudimentary way of doing a sim. run, with a million variates queued up: len = 15 box = 2^len - 1 cnt = 1 Catch[Scan[(If[(box = BitXor[box, #]) == 0, Throw[cnt], cnt++]) &, 2^RandomInteger[len - 1, 1000000]]] Folding it instead of scanning may well be faster. If you want to solve directly, say for the case of a box of six: d = DiscreteMarkovProcess[1, { {0, 1, 0, 0, 0, 0, 0}, {1/6, 0, 5/6, 0, 0, 0, 0}, {0, 2/6, 0, 4/6, 0, 0, 0}, {0, 0, 3/6, 0, 3/6, 0, 0}, {0, 0, 0, 4/6, 0, 2/6, 0}, {0, 0, 0, 0, 5/6, 0, 1/6}, {0, 0, 0, 0, 0, 0, 1} }]; Mean[FirstPassageTimeDistribution[d, 7]] (* 416/5 *) For the general solution: size[n_] := Module[{mp, sa}, sa = SparseArray[{Band[{2, 1}] -> Range[1/n, (n - 1)/n, 1/n], Band[{2, 3}] -> Range[(n - 1)/n, 1/n, -1/n], {1, 2} -> 1, {n + 1, n + 1} -> 1}, {n + 1, n + 1}]; mp = DiscreteMarkovProcess[1, sa]; Mean[FirstPassageTimeDistribution[mp, n + 1]]]; size[20] (* 3234139734016/2909907 *) Which takes only a fraction of a second even on the netbook I'm on right now.