Barnsley Fern via R

Here's the deal. Start with the initial point \( p_0=\begin{matrix} 0 \\ 0 \end{matrix} \). Then, for each iteration, the point get transformed in the following fashion. With 1% probability, \[ p_{n+1}=\begin{matrix} 0 & 0 \\ 0 & 0.16 \end{matrix} p_n \] With 85% probability, \[ p_{n+1}=\begin{matrix} 0.85 & 0.04 \\ -0.04 & 0.85 \end{matrix} p_n + \begin{matrix} 0 \\ 0.16 \end{matrix} \] With 7% probability, \[ p_{n+1}=\begin{matrix} 0.2 & -0.26 \\ 0.23 & 0.22 \end{matrix} p_n +\begin{matrix} 0 \\ 0.16 \end{matrix} \] and with the rest 7%, \[ p_{n+1}=\begin{matrix} -0.15 & 0.28 \\ 0.26 & 0.24 \end{matrix} p_n +\begin{matrix} 0 \\ 0.44 \end{matrix} \]

iter = 10000
p = runif(iter)
coord = c(0, 0)
df = data.frame(t(coord))
for (i in 1:iter) {
    if (p[i] <= 0.01) {
        m = matrix(c(0, 0, 0, 0.16), nrow = 2, ncol = 2)
        const = c(0, 0)
    } else if (p[i] > 0.01 && p[i] <= 0.86) {
        m = matrix(c(0.85, -0.04, 0.04, 0.85), nrow = 2, ncol = 2)
        const = c(0, 1.6)
    } else if (p[i] > 0.86 && p[i] <= 0.93) {
        m = matrix(c(0.2, 0.23, -0.26, 0.22), nrow = 2, ncol = 2)
        const = c(0, 1.6)
    } else {
        m = matrix(c(-0.15, 0.26, 0.28, 0.24), nrow = 2, ncol = 2)
        const = c(0, 0.44)
    }
    coord = m %*% coord + const
    df = rbind(df, c(coord))
}

plot(x = df[, 2], y = df[, 1], plt = c(0, 10, -5, 5), cex = 0.1, asp = 1)