Algorithm to estimate the probability that a "0 and 1 matrix" fills up following the bootstrap percolation rules Presentation of the model: we consider the regular lattice created from $\mathbb{Z}^2$ (It's no more, no less a square lattice).


At $t=0$, each site is said "active" independently with a probability $p$, "inactive" otherwise.
At $t$, if a site is inactive, it becomes active if at least $2$ of its neighbours are active.


We call neighbour of $x$ each element of $V(x)=\{y\in\mathbb{Z}^2, ||x-y||=1\}$ (or, if I have a matrix M in python, a neighbour of $M[i][j]$ is $M[i\pm1][j]$ or $M[i][j\pm 1]$).



My problem: I consider a $n\times n$ matrix filled with 0s and 1s, (0 = inactive ; 1 = active), and I want to estimate the probability $p_c$ that a matrix, following the rules above, fills up (at $t\rightarrow \infty$, but this happen really fast in reality). The probability $p_c$ depend on $n$. I want to test it for $n=2$ to $n=100$ at least.

What I already have: I have a program, matrix_creation(n,p) that creates a $n\times n$  matrix with the probability $p$ for each coefficient to be 1, and another, test(M), to test if "M" fills up.
Here is a reasonably fast approach to estimate $p_c$. We will perform $N$ different trials, count the number of successes $M$, and estimate the result as $M/N \pm 1/\sqrt{N}$ (better estimates are possible if $p_c$ is close to 0 or to 1).

For each trial, you initialize an $N \times N$ matrix with the initial filling. You then go over the array, calculating the number of active neighbors of each non-filled point in the matrix, and whenever that number reaches 2, you put it in a stack. You then repeatedly pop elements from the stack, updating the matrix of active neighbors, and whenever that number reaches 2 for some non-filled point, you add it to the stack. When the stack is empty, the process is complete, and you can check whether all points have been activated.

There are many optimizations possible: for example, you can use an $(N+2)\times(N+2)$ matrix to help deal with border cases, you can keep track of the number of active points, and so on.

Using this naive approach, I managed to perform $10^6$ trials for $n = 100$, $p = 0.0535$ in about two minutes on my PC. The answer was $M/N \approx 0.5008$, so the critical probability is pretty close to $0.0535$.