I have figured out how to do the sierpinski procedure for a triangle but not for a square, which i need.

serpinski:=proc(n::posint)
PLOT(CURVES([serp([0,0],[1/2.,sqrt(3.0)/2.0],[1,0],n)]),AXESSTYLE(NONE),SCALING(CONSTRAINED));
end:
 

serp:=proc(p,q,r,s,n)
local m1,m2,m3,m4;
if n=0 then return(p,q,r,s,p)
end if;
m1:=1/2.*(p+q);
m2:=1/2.*(q+r);
m3:=1/2.*(r+s);
m4:=1/2.*(s+p);
serp(p,m1,m4,n-1),m1,serp(m1,q,m2,n-1),m3,serp(m3,m2,r,n-1),m4,serp(p
end:
 

this is what i have for the triangle...