Barnsley fern

BMat=rbind(c(0.849,0.037,-0.037,0.849,0.075,0.183),c(0.197,-0.226,0.226,0.197,0.400,0.049),c(-0.150,0.283,0.260,0.237,0.575,-0.084),c(0,0,0,0.160,0.500,0))

# Initial conditions:
x=0
y=0
X11()
plot(0,0,xlim=c(0,1),ylim=c(0,1),col="white")
COLOR=c("green","red","blue","yellow")

for(j in 1:10000)
{
i=sample(1:4,1) # ,prob=c(0.25,0.25,0.25,0.25)
x3=x
x=BMat[i,1]*x+BMat[i,2]*y+BMat[i,5]
y=BMat[i,3]*x3+BMat[i,4]*y+BMat[i,6]
points(x,y,pch=".",cex=1, col=COLOR[i])
}

#--------------------------------------#
# Alternative code by Stefano Guazzetti
#--------------------------------------#

# The probabilities used here give more dense fern.

niter<-150000
kind<-sample(1:4, niter, repl=T, prob=c(.01, .07, .07, .85))

x<-numeric(niter+1)
y<-numeric(niter+1)
x[1]<-0
y[1]<-0
for (i in 1:niter) {
x[i+1]<- 0*(kind[i]==1)+(0.2*x[i]-0.26*y[i])*(kind[i]==2)+
(-0.15*x[i]+0.28*y[i])*(kind[i]==3) +
(0.85*x[i]+0.04*y[i])*(kind[i]==4)

y[i+1]<- 0.16*y[i]*(kind[i]==1)+(0.23*x[i]+0.22*y[i]+1.6)*(kind[i]==2)+
(0.26*x[i]+0.24*y[i]+.44)*(kind[i]==3) +
(-0.04*x[i]+0.85*y[i]+1.6)*(kind[i]==4)
}

par(mar=c(0.1,0.1,0.1,0.1))
plot(x, y, pch=17, cex=.3, col="darkgreen", axes=F, ann=F)