restart;
with(IterativeMaps):
form := z - z^2;

                                     2
                           form := -z  + z

z := subs(z=zr+zi*I, form):
fzi := evalc( Im( z ) );

                        fzi := -2 zi zr + zi

fzr := evalc( Re( subs(zi=zitemp, z) ) );

                                   2     2
                      fzr := zitemp  - zr  + zr

julia_img := CodeTools:-Usage(
    Escape( [zi  , zr  , zitemp , zrsqr , zisqr],
            [fzi , fzr , zi     , zr^2  , zi^2 ],
            [y   , x   , y      , x^2   , y^2  ],
            zrsqr+zisqr > 250,
            -1.0, 2.0, -1.5, 1.5,
            iterations = 30,
            height=600, width=600 ) ):

    memory used=20.41MiB, alloc change=46.80MiB, cpu time=656.00ms,
    real time=285.00ms, gc time=0ns

with(ImageTools):
FitIntensity(julia_img, inplace):
Embed(julia_img);

restart;
form := z - z^2;

                          2
                        -z  + z

update := op(evalc([Re,Im](subs(z=re+im*I, form))));

                           2     2
               update := im  - re  + re, -2 im re + im

JuliaSet := subs(__dummy=update,
  proc(a, b)
    local re, im, resq, imsq, m;
    (re, im) := (a, b):
    resq := re^2:
    imsq := im^2;
    for m to 30 while resq+imsq < 250 do
      (re,im) := __dummy;
      resq := re^2;
      imsq := im^2;
    end do;
    return m;
  end proc):

CodeTools:-Usage(
  plots:-densityplot(JuliaSet, -1.0..2.0, -1.5..1.5,
                     colorstyle=HUE, grid=[600, 600],
                     style=patchnogrid, axes=none,
                     size=[600,600]) );

    memory used=4.93GiB, alloc change=8.24MiB, cpu time=25.88s,
    real time=25.90s, gc time=5.32s