I made this with SAGE, an open-source math package, for entertainment, as an exercise. Permalink here; original source code (not this version) follows:

#*****************************************************************************
#       Copyright (C) 2008 Dean Moore
#
#
#  Distributed under the terms of the GNU General Public License (GPL)
#                  http://www.gnu.org/licenses/
#*****************************************************************************
# Animates the construction of the Butterfly Curve.
# Written by Dean Moore, March 30 -- April 01, 2008

def animate_curve((g,f), a, b, step, x_min, x_max, y_min, y_max, fig_size = 5):
    v = []
    def c(i):                                 # For output
       return ((i-a)/(b-a), 0, 1-(i-a)/(b-a)) # colors.

    for i in srange(a, b + step, step):
        if i >= b:
            curve_color=(0,0,1) # Final image is blue.
        else:
            curve_color = c(i)

        p = parametric_plot((g,f), 0, i, rgbcolor=curve_color, plot_points=50*(i-a+1)/step, thickness = 0.25)
        if i < b:
            p = p + points((g(i),f(i)),rgbcolor=(0,0,0), pointsize=10) # Terminal point on all but last image.
        v.append(p)
    a = animate(v, xmin=x_min, xmax = x_max, ymin = y_min, ymax = y_max, figsize=[fig_size,fig_size])
    return a

figure_size                            = 4
x_min = -4;x_max = 4;y_min = -2; y_max = 5
number_of_images_in_final_frozen_image = 12
delay_between_images                   = 0 

def g(x):
    return sin(x)*(e^(cos(x)) - 2*cos(4*x) - (sin(x/12))^5) # Parametric 
def f(x):                                                   # equations of the
    return cos(x)*(e^(cos(x)) - 2*cos(4*x) - (sin(x/12))^5) # Butterfly Curve. 

step = 0.5
b = animate_curve((g,f),0, 24*pi, step, x_min, x_max, y_min, y_max, figure_size)

d = animate([b[-1] for i in srange(0, number_of_images_in_final_frozen_image)], # Last image
             xmin=x_min, ymin=y_min,                                            # is "held" a
             xmax=x_max, ymax=y_max,                                            # time.
             figsize=[figure_size,figure_size])