This isn't the fastest mandelbrot set algorithm, but by far the coolest I've yet written. Creating a discrete 2d space containing complex floats and mapping the function on it. Unfortunately I did not find a builtin way to obtain that space from ranges. The complex numbers seem to create problems. But I'm only starting :)

using PyPlot

function f(x0)
    nmax = 50
    xmax = 4
    c = x0
    x = x0
    n = 0

    while abs(x) < 4 && n < nmax
        x = x * x + c
        n += 1
    end

    n
end

function cspace(ca, cb, cn)
    nx = real(cn)
    ny = imag(cn)
    dx = (real(cb)-real(ca)) / nx
    dy = (imag(cb)-imag(ca)) / ny

    [complex(real(ca) + i*dx, imag(ca) + j*dy) for i in 0:nx, j in 0:ny]
end

matshow(f.(cspace(-3-1im, 1+1im, 200+200im)))