Here is an attempt with MetaPost, for whom it may interest. The recursive macro (square_fractal
) at the basis of this program is heavily inspired by this answer to a closely related subject.
vardef square_fractal(expr A, B, n) =
save P; pair P[]; P0 = A; P1 = B;
for i = 1 upto 2:
P[i+1] = P[i-1] rotatedaround (P[i], -90);
endfor;
if n = 0: fill P0 for i = 1 upto 3: -- P[i] endfor -- cycle;
else:
save Q; pair Q[];
for i = 0, 2:
Q[i] = 1/4[P[i],P[i+1]]; Q[i+1] = 3/4[P[i],P[i+1]];
square_fractal(P[i], Q[i], n-1);
square_fractal(Q[i+1], P[i+1], n-1);
endfor;
square_fractal(P0 rotatedaround (Q0, -90), P1 rotatedaround (Q1, 90), n-1); fi
enddef;
beginfig(1);
for n = 0 upto 4:
draw image(square_fractal(origin, (4cm, 0), n)) shifted (n*4.5cm, 0);
endfor;
endfig;
end.

Starting from order 0 (the full square), MetaPost manages an output up to order 6 on my machine. Interestingly enough, order 7 is reached if the previous code is included in a LuaLaTeX program. I don't know the reason why.
Edit Still within LuaLaTeX, and after using floating point numerics (\mplibnumbersystem{double}
added right after \usepackage{luamplib}
) instead of the defaults fixed point numerics , MetaPost manages to produce the figure at order 9 after 20 minutes. But it nearly freezes my very old laptop (a MacBook Pro from 2008), so I don't dare to go further on it. Maybe I will try it again on a more recent and more powerful computer.
\RequirePackage{luatex85}
\documentclass[border=5mm]{standalone}
\usepackage{luamplib}
\mplibnumbersystem{double}
\begin{document}
\begin{mplibcode}
vardef square_fractal(expr A, B, n) =
save P; pair P[]; P0 = A; P1 = B;
for i = 1 upto 2:
P[i+1] = P[i-1] rotatedaround (P[i], -90);
endfor;
if n = 0: fill P0 for i = 1 upto 3: -- P[i] endfor -- cycle;
else:
save Q; pair Q[];
for i = 0, 2:
Q[i] = 1/4[P[i],P[i+1]]; Q[i+1] = 3/4[P[i],P[i+1]];
square_fractal(P[i], Q[i], n-1);
square_fractal(Q[i+1], P[i+1], n-1);
endfor;
square_fractal(P0 rotatedaround (Q0, -90), P1 rotatedaround (Q1, 90), n-1); fi
enddef;
beginfig(1);
square_fractal(origin, (12cm, 0), 9);
endfig;
\end{mplibcode}
\end{document}
The figure below is the one of order 8. I couldn't manage to produce a PNG version of order 9 because of the near-freezing of my laptop.
