Code

The main commands are GraphicsIFS[system,starting picture,number of iterations] and RandomIFS[system,number of points] together with the "try a system" commands which are:
	Curlicue[]
	SierpinskiTriangle[]
	SierpinskiCarpet[]
	CantorSquare[]
Also one can generate a picture of your own system with AskASystem[number of maps]

Off[General::spell]
Off[General::spell1]

<< LinearAlgebra`MatrixManipulation`

Clear[NMaps, Maps, Dim]

    
NMaps[IFS[nmaps_Integer, As_List, bs_List]] := 
        nmaps;

    
Maps[IFS[nmaps_Integer, As_List, bs_List], m_Integer, pos_List] := 
      As[[m]] . pos + 
        bs[[m]]

    
Dim[IFS[nmaps_Integer, As_List, bs_List]]  := 
      Length[bs[[1]]]

Clear[DoRandomMap, RandomIFS, Converter, 
      RandomWithDistribution]

    
Converter[probs_List]  :=  
      FoldList[Plus, 0, probs]
RandomWithDistribution[dist_List]  :=  
      Module[{rand}, rand = Random[]; 
        Length[Select[dist, (rand > #  &)]]]
DoRandomMap[sys_IFS, pos_List]  := 
      Maps[sys, Random[Integer, {1, NMaps[sys]}], 
        pos]
DoRandomMap[sys_IFS, dist_List, pos_List]  := 
      Maps[sys, RandomWithDistribution[dist], pos]

    
Options[RandomIFS]  =  {DropFirst -> 1, Probs -> {}, 
          DropInterval -> {}};
RandomIFS[sys_IFS, steps_Integer, opts___]  :=  
      Module[{probs, dropfirst, repeater, dropcoords, dropper}, 
        dropfirst = DropFirst  /.  {opts}  /.  Options[RandomIFS]; 
        probs = Probs /. {opts} /. Options[RandomIFS]; 
        dropcoords  =  DropInterval  /.  {opts}  /.  
            Options[RandomIFS]; 
        If[probs == {}, repeater = (DoRandomMap[sys, #]  &), 
          repeater = (DoRandomMap[sys, Converter[probs], #] &)]; 
        If[dropcoords == {}, 
          ShowStraight[
            Graphics[
              Point  /@  
                Drop[NestList[repeater, Table[Random[], {2}], steps], 
                  dropfirst]]], dropper = (Drop[#, dropcoords]  &); 
          ShowStraight[
            Graphics[
              Point  /@  
                Drop[dropper  /@  
                    NestList[repeater, Table[Random[], {Dim[sys]}], steps], 
                  dropfirst]]]]]

Clear[ApplyIFS, ApplyMap, DropFrom, GraphicsIFS, 
        ShowStraight];

    
DropFrom[Line[l_List], drop_]  :=  
      Line[(Drop[#, drop] &) /@ l]
DropFrom[Polygon[l_List], drop_] := 
      Polygon[(Drop[#, drop] &) /@ l]
DropFrom[Point[l_List], drop_] := 
      Point[(Drop[#, drop] &) /@ l]

    
ApplyMap[sys_IFS, m_Integer, Line[poss_List]]  :=  
      Line[(Maps[sys, m, #] &) /@ poss]
ApplyMap[sys_IFS, m_Integer, Polygon[poss_List]] := 
      Polygon[(Maps[sys, m, #] &) /@ poss]
ApplyMap[sys_IFS, m_Integer, Point[pos_List]] := 
      Point[Maps[sys, m, pos]]
ApplyMap[sys_IFS, m_Integer, 
        l_List] := (ApplyMap[sys, m, #]  &) /@  
        l

    
ApplyIFS[sys_IFS, l_List]  :=  
      Join  @@  (ApplyIFS[sys, #]  &)  /@  
          l
ApplyIFS[sys_IFS, item_]  :=  
      Table[ApplyMap[sys, m, item], {m, NMaps[sys]}]

    
Options[
          GraphicsIFS] = {DropInterval -> {}};
GraphicsIFS[sys_IFS, init_, steps_Integer, opts___]  :=  
      Module[{dropcoords}, 
        dropcoords  =  DropInterval  /.  {opts}  /.  
            Options[RandomIFS]; 
        If[dropcoords == {}, 
          Graphics[Nest[(ApplyIFS[sys, #]  &), init, steps]]  // 
            ShowStraight, dropper  =  (DropFrom[#, dropcoords] &); 
          Graphics[
              dropper  /@  
                Nest[(ApplyIFS[sys, #]  &), init, steps]]  // 
            ShowStraight]]

    
Graphics3DIFS[sys_IFS, init_, steps_Integer]  :=  
      Graphics3D[Nest[(ApplyIFS[sys, #]  &), init, steps]]  //  
        ShowStraight

    
ShowStraight[args__]  :=  
      Show[args, AspectRatio -> 1, PlotRange -> All]

sq  =  
        Line[{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 
              0}}];
cross  =  {Line[{{0.5, 0}, {0.5, 1}}], 
          Line[{{0, 0.5}, {1, 0.5}}]};
nhappy  =  20;
happy  =  {Line[
            Table[{0.5, 0.5} + 
                0.5 {Cos[2  Pi i/nhappy], Sin[2  Pi i/nhappy]}, {i, 0, 
                nhappy}]], 
          Line[Table[{0.5, 0.5} + 
                0.25 {Cos[2  Pi i/nhappy], Sin[2  Pi i/nhappy]}, {i, 
                nhappy/2, nhappy}]], 
          Line[Table[{0.3, 0.7} + 
                0.1 {Cos[2  Pi i/nhappy], Sin[2  Pi i/nhappy]}, {i, 0, 
                nhappy}]], 
          Line[Table[{0.7, 0.7} + 
                0.1 {Cos[2  Pi i/nhappy], Sin[2  Pi i/nhappy]}, {i, 0, 
                nhappy}]]};
block  =  
        Polygon[{{0, 0}, {1, 0}, {1, 1}, {0, 1}, {0, 
              0}}];
horiz = Line[{{0, 0}, {1, 0}}];
spike = 
        Line[{{0, 0}, {0.5, 1}, {1, 0}}];
diag = Line[{{0, 0}, {1, 1}}];
flatline = Line[{{0, 0, 0}, {1, 1, 0}}];
bentline = 
        Line[{{0, 0, 0}, {0.5, 0.5, 1}, {1, 1, 
              0}}];
diagline = Line[{{0, 0, 0}, {1, 1, 1}}];

Clear[Interpolating, Threed];

    
Interpolating[x_List, vscales_List]  := 
      Module[{a, e, c, f, n, i}, n = Length[x]; 
        a = Table[
              x[[i + 1, 1]] - 
                x[[i, 1]], {i, 
                n - 1}]/(x[[n, 1]] - 
                x[[1, 1]]); 
        e = Table[
              x[[n, 1]] 
                  x[[i - 1, 1]] - 
                x[[1, 1]] 
                  x[[i, 1]], {i, 2, 
                n}]/(x[[n, 1]] - 
                x[[1, 1]]); 
        c = Table[
              x[[i + 1, 2]] - 
                x[[i, 2]] - 
                vscales[[
                    i]] (x[[n, 
                        2]] - 
                      x[[1, 2]]), {i, 
                n - 1}]/(x[[n, 1]] - 
                x[[1, 1]]); 
        f = Table[
              x[[n, 1]] 
                  x[[i - 1, 2]] - 
                x[[1, 1]] 
                  x[[i, 2]] - 
                vscales[[
                    i - 1]] (x[[n, 
                          1]] 
                        x[[1, 2]] - 
                      x[[1, 1]] 
                        x[[n, 
                          2]]), {i, 2, 
                n}]/(x[[n, 1]] - 
                x[[1, 1]]); 
        Print[a, c, e, f]; 
        IFS[n - 1, 
          Table[{{a[[i]], 
                0}, {c[[i]], 
                vscales[[i]]}}, {i, 
              n - 1}], 
          Table[{e[[i]], 
              f[[i]]}, {i, 
              n - 1}]]]

    
Threed[x_List, H_List, vmaps_List]  := 
      Module[{a, c, k, e, f, g, n, xrange, Frange, Hrange}, n = Length[x]; 
        xrange = 
          x[[n, 1]] - 
            x[[1, 1]]; 
        Frange = 
          x[[n, 2]] - 
            x[[1, 2]]; 
        Hrange = 
          H[[n]] - 
            H[[1]]; 
        a = Table[
              x[[i + 1, 1]] - 
                x[[i, 1]], {i, n - 1}]/
            xrange; c = 
          Table[x[[i + 1, 2]] - 
                x[[i, 2]] - 
                vmaps[[i, 1, 1]]*
                  Frange - 
                vmaps[[i, 1, 2]]*
                  Hrange, {i, n - 1}]/xrange; 
        k = Table[
              H[[i + 1]] - 
                H[[i]] - 
                vmaps[[i, 2, 1]]*
                  Frange - 
                vmaps[[i, 2, 2]]*
                  Hrange, {i, n - 1}]/xrange;  
        e = Table[
              x[[n, 1]] 
                  x[[i, 1]] - 
                x[[1, 1]] 
                  x[[i + 1, 1]], {i, 
                n - 1}]/xrange; 
        f = Table[
            x[[i, 2]] - 
              vmaps[[i, 1, 1]] 
                x[[1, 2]] - 
              vmaps[[i, 1, 2]] 
                H[[1]] - 
              c[[i]] 
                x[[1, 1]], {i, n - 1}]; 
        g = Table[
            H[[i]] - 
              vmaps[[i, 2, 1]] 
                x[[1, 2]] - 
              vmaps[[i, 2, 2]] 
                H[[1]] - 
              k[[i]] 
                x[[1, 1]], {i, n - 1}]; 
        IFS[n - 1, 
          Table[{{a[[i]], 0, 
                0}, {c[[i]], 
                vmaps[[i, 1, 1]], 
                vmaps[[i, 1, 
                  2]]}, {k[[i]], vmaps[[i, 2, 1]], 
                vmaps[[i, 2, 2]]}}, {i, 
              n - 1}], 
          Table[{e[[i]], 
              f[[i]], 
              g[[i]]}, {i, 
              n - 1}]]]

Clear[AskASystem, CantorSquare, SierpinskiCarpet, SierpinskiTriangle, 
      Curlicue, 
      ExtractConstant]

\(\(strings\  = \ {"first", "second", "third", "fourth", "fifth", "sixth", "seventh", "eighth", "ninth", "tenth"};\)
\)
ExtractConstant[poly_]  :=  
      Module[{coeffs}, coeffs  =  CoefficientList[poly, {x, y}]; 
        If[Length[coeffs] == 0, Return[0], 
          If[Length[
                coeffs[[1]]] == 0, 
            0, coeffs[[1, 
              1]]]]]

\(AskASystem[n_]\  := \ 
      Module[{inputs, coeffs, A, b, sys}, inputs\  = \ Table[0, {n}, {2}]; 
        coeffs = Table[0, {2}]; A = Table[0, {n}]; b = Table[0, {n}]; 
        For[i = 1, i ≤ n, \(i++\), 
          For[j = 1, j ≤ 2, \(j++\), 
            inputs[[i, j]] = 
              Input["What is the " <> 
                  strings[[
                    j]] <> " component of the " <> 
                  strings[[
                    i]] <> " map?"]]; 
A[[
              i]] = {{Coefficient[
                  inputs[[i, 1]], x], 
                Coefficient[
                  inputs[[i, 1]], 
                  y]}, {Coefficient[
                  inputs[[i, 2]], x], 
                Coefficient[
                  inputs[[i, 2]], y]}}; 
          b[[i]] = {ExtractConstant[
                inputs[[i, 1]]], 
              ExtractConstant[
                inputs[[i, 2]]]}]; 
        sys = IFS[n, A, b]; 
        For[i = 1, i ≤ 5, \(i++\), GraphicsIFS[sys, block, i]; 
          Pause[3]]; 
        GraphicsIFS[sys, block, 
          6]]
\)
SierpinskiTriangle[]  :=  AskASystem[3]
CantorSquare[] := AskASystem[4]
SierpinskiCarpet[] := AskASystem[8]
Curlicue[] := AskASystem[3]

ExtractConstant[0]

Curlicue[]

Systems

The systems included are:
	The Sierpinksi triangle	sierp
    	The Cantor dust		cantor
    	An L-shaped carpet		carpet
    	A design with lots of hooks	hooks
    	A mystery fractal		mystery
    	
 The intial shapes are:
 	A block			block
 	A square			sq
 	A smiley face			happy
 	A cross				cross

Clear[Rot];

    
Rot[θ_]  :=  {{Cos[θ], -Sin[θ]}, {Sin[θ], Cos[θ]}}

A[1]  =  {{1/2, 0}, {0, 1/2}};
b[1]  =  {0, 0};
A[2]  =  A[1];
b[2]  =  {1/2, 0};
A[3]  =  A[1];
b[3]  =  {1/4, 1/2};

    
sierp  =  
      IFS[3, Table[A[i], {i, 3}], 
        Table[b[i], {i, 3}]]

A[1]  =  {{1/3, 0}, {0, 1/3}};
b[1]  =  {0, 0};
A[2]  =  A[1];
b[2]  =  {2/3, 0};
A[3]  =  A[1];
b[3]  =  {2/3, 2/3};
A[4]  =  A[1];
b[4]  =  {0, 2/3};

    
cantor  =  
      IFS[4, Table[A[i], {i, 4}], Table[b[i], {i, 4}]]

A[
          1]  =  {{0.849, 0.037}, {-0.037, 
            0.849}};
b[1]  =  {0.075, 0.183};
A[
          2]  =  {{0.197, -0.226}, {0.226, 
            0.197}};
b[2]  =  {0.4, 0.0490};
A[
          3]  =  {{-0.15, 0.283}, {0.260, 
            0.237}};
b[3]  =  {0.575, -0.084};
A[4]  =  {{0.0, 0.0}, {0.0, 0.16}};
b[4]  =  {0.5, 0};

mystery  =  
      IFS[4, Table[A[i], {i, 4}], 
        Table[b[i], {i, 4}]]

A[1]  =  {{0.5, 0}, {0, 0.5}};
b[1]  =  {0, 0};
A[2]  =  {{0.5, 0}, {0, 0.5}};
b[2]  =  {0.5, 0};
A[3]  =  {{0, -0.5}, {0.5, 0}};
b[3]  =  {1, 0.5};

    
hooks  =  
      IFS[3, Table[A[i], {i, 3}], 
        Table[b[i], {i, 3}]]

A[1]  =  (0.5) Rot[Pi/2] ;
b[1] = {0.5, 0};
A[2]  =  (0.5) Rot[0];
b[2]  =  {0, 0.5};
A[3]  =  (0.5) Rot[-Pi/2];
b[3] = {0.5, 1};
μ = 3/8;
A[4]  =  μ Rot[0]  //  N;
b[4] =  {0.5, 0.5} - A[4] . {0.5, 0.5};

    
carpet  =  
      IFS[4, Table[A[i], {i, 4}], Table[b[i], {i, 4}]]

A[1]  =  0.5 {{0, 1}, {1, 0}};
b[1] = {0, 0};
A[2] = 0.5  IdentityMatrix[2];
b[2] = {0, 0.5};
A[3] = A[2];
b[3] = {0.5, 0.5};
A[4] = -A[1];
b[4] = {1, 0.5};

    
spacey = 
      IFS[4, Table[A[i], {i, 4}], 
        Table[b[i], {i, 4}]]
hiddenspacey  =  
      IFS[4, Table[
          BlockMatrix[{{{{0.25}}, {{0, 0}}}, {{{0}, {0}}, A[i]}}], {i, 4}], 
        Table[Join[{(i - 1.0)/4}, b[i]], {i, 4}]]

On[General::spell]
On[General::spell1]

Examples

Curlicue[]

GraphicsIFS[sierp, block, 1]

RandomIFS[sierp, 100]

RandomIFS[fern, 1000, Probs -> {0.7, 0.1, 0.1, 0.1}, 
      DropFirst -> 20]