The following Maxima code was sent to the Maxima mailing list by Jaime Villate <villate@fe.up.pt>. It takes a little while to run (at least on my machine), so give it some time. Enjoy! a1: matrix([0.85,0],[0,0.85])$ a2: matrix([-0.2,0.26],[-0.23,-0.22])$ a3: matrix([0.15,-0.28],[-0.26,-0.24])$ a4: matrix([0,0],[0,0.16])$ q1: [0,1.6]$ q2: [0,0]$ m1: matrix([1/2,sqrt(3)/6],[-sqrt(3)/6,1/2])$ m2: matrix([1/2,-sqrt(3)/6],[sqrt(3)/6,1/2])$ l: [[-0.4,11],[0,12.1],[0.4,11],[-0.4,11]]$ /* Koch's fractal */ for j thru 3 do (l2: [l[1]], for i thru length(l)-1 do ( p1: l[i], p2: l[i+1], l2: endcons(float((2*p1+p2)/3), l2), l2: endcons(float(first(transpose(m1.p1+m2.p2))), l2), l2: endcons(float((p1+2*p2)/3), l2), l2: endcons(p2,l2)), l: l2)$ /* Barnsley's Iterative Function System */ ifs(prob, mat, point, p0, n) := block ([p:p0, s, r, m:last(prob), data:[p0]], for i thru n do (r: random(m) + 1, s: 0, for j while (r-prob[j]) > 0 do s:j, p: float(mat[s+1].p + point[s+1]), data: cons(first(transpose(p)), data)), data)$ l1: ifs([85,92,99,100], [a1,a2,a3,a4], [q1,q1,q1,q2], [0,0], 30000)$ plot2d([[discrete,l],[discrete,l1]],[style,lines,dots],[legend,false],[color,blue,green],[box,false],[x,-8,8]);