auto <- function(n = 500) {
    a <- matrix(0, n, n)
    a[1, 1] <- 1
    for (i in 2:n) a[i, 2:n] <- a[i - 1, 1:(n - 1)] + a[i - 1, 2:n] == 1
    plot(expand.grid(1:n, 1:n), col = apply(a, 2, rev), asp = 1, pch = ".", xaxt = "n", yaxt="n", bty="n",ylab="",xlab="")
}

for(i in seq(3, 500, length = 10))auto(i)