Algorithm to estimate the probability that a “0 and 1 matrix” fills up following the bootstrap percolation rules

Question

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.