code

> basic := proc(p1,p2,p3)
local dx,dy,p4,p5,p6,p7,p8,p9;
p4 := .5*(p1+p2);
p9 := .5*(p2+p3);
p5 := p4+(p2-p9);
p6 := p2+(p2-p9);
p7 := p2+(p4-p1);
p8 := p9+(p4-p1);
p4,p5,p6, p2 ,p7,p8,p9, p3 ;
end:

>

> peano := proc(fl)
local i,cur;
cur := [fl[1] ] ;
for i from 1 by 2 to nops(fl)-2 do
cur := [op(cur),basic( fl[i],fl[i+1],fl[i+2] )]
od;

> end:

> sqr := plottools[polygon]([[0,0],[2,0],[2,2],[0,2]],color=yellow):

> hilbert := proc(n,clr)
local fl,i,ti;

> fl := [[0,0],[1,1],[2,0]]:

> for i from 1 to n-1 do fl := peano(fl) od:
ti := cat(n,`th stage of Peano's curve`);
plots[display]([sqr,plot(fl,style=LINE,axes=NONE,color=clr,
thickness=3,title=ti,
scaling=CONSTRAINED)]); end:

>

>

> clrs := [blue,red,black,tan,magenta]:
plots[display]([seq(hilbert(i,clrs[i]),i=1..5)],insequence=true,scaling=constrained);