P2B2:=proc(z1,z2,z3,d)    a[0]:=z1: b[0]:=z2: c[0]:=z3:
>     for h from 1 to d do
>     a[h]:=a[h-1]: b[h]:=(1/3-1/sqrt(3))*a[h-1]+(1/3)*b[h-1]-(2/3)*c[h-1]: c[h]:=(1/3-1/sqrt(3))*a[h-1]-(2/3)*b[h-1]+(1/3)*c[h-1]:
>     unassign('a[h-1]','b[h-1]','c[h-1]'):
>     od:
>     for i from 0 to d-1 do
>     theta[2*i]:=3+sqrt(3): theta[2*i+1]:=3-sqrt(3):
>     od:
>     X:=array(1..(1+2*d)): Y:=array(2..(1+2*d)):
>     X:=array(1..(2*d+1)): Y:=array(2..(2*d+1)):
>     X[1]:=array(1..1): X[1][1]:=array(1..3): X[1][1][1]:=array(1..3,[0,0,a[d]]): X[1][1][2]:=array(1..3,[4^d,0,b[d]]): X[1][1][3]:=array(1..3,[(4^d)/2,4^d*(sqrt(3))/2,c[d]]):
>     for k from 2 to (2*d+1) do
>     X[k]:=array(1..(3^(k-1))): Y[k]:=array(1..(3^(k-2))):
>     for l from 1 to (3^(k-2)) do
>     X[k][3*l-2]:=array(1..3): X[k][3*l-1]:=array(1..3): X[k][3*l]:=array(1..3):        Y[k][l]:=array(1..6):
>     Y[k][l][1]:=array(1..3,[X[k-1][l][1][1],X[k-1][l][1][2],X[k-1][l][1][3]]): Y[k][l][3]:=array(1..3,[X[k-1][l][2][1],X[k-1][l][2][2],X[k-1][l][2][3]]): Y[k][l][5]:=array(1..3,[X[k-1][l][3][1],X[k-1][l][3][2],X[k-1][l][3][3]]): Y[k][l][2]:=array(1..3,[(X[k-1][l][1][1]+X[k-1][l][2][1])/2,X[k-1][l][1][2],((4-theta[2*d+1-k])*(X[k-1][l][1][3]+X[k-1][l][2][3])+2*X[k-1][l][3][3])/((2-theta[2*d+1-k])*(5-theta[2*d+1-k]))]): Y[k][l][4]:=array(1..3,[(X[k-1][l][2][1]+X[k-1][l][3][1])/2,(X[k-1][l][2][2]+X[k-1][l][3][2])/2,((4-theta[2*d+1-k])*(X[k-1][l][2][3]+X[k-1][l][3][3])+2*X[k-1][l][1][3])/((2-theta[2*d+1-k])*(5-theta[2*d+1-k]))]): Y[k][l][6]:=array(1..3,[(X[k-1][l][3][1]+X[k-1][l][1][1])/2,(X[k-1][l][3][2]+X[k-1][l][1][2])/2,((4-theta[2*d+1-k])*(X[k-1][l][3][3]+X[k-1][l][1][3])+2*X[k-1][l][2][3])/((2-theta[2*d+1-k])*(5-theta[2*d+1-k]))]):
>     X[k][3*l-2]:=array(1..3): X[k][3*l-2][1]:=array(1..3,[Y[k][l][1][1],Y[k][l][1][2],Y[k][l][1][3]]): X[k][3*l-2][2]:=array(1..3,[Y[k][l][2][1],Y[k][l][2][2],Y[k][l][2][3]]): X[k][3*l-2][3]:=array(1..3,[Y[k][l][6][1],Y[k][l][6][2],Y[k][l][6][3]]):     X[k][3*l-1]:=array(1..3): X[k][3*l-1][1]:=array(1..3,[Y[k][l][2][1],Y[k][l][2][2],Y[k][l][2][3]]): X[k][3*l-1][2]:=array(1..3,[Y[k][l][3][1],Y[k][l][3][2],Y[k][l][3][3]]): X[k][3*l-1][3]:=array(1..3,[Y[k][l][4][1],Y[k][l][4][2],Y[k][l][4][3]]): X[k][3*l]:=array(1..3): X[k][3*l][1]:=array(1..3,[Y[k][l][6][1],Y[k][l][6][2],Y[k][l][6][3]]): X[k][3*l][2]:=array(1..3,[Y[k][l][4][1],Y[k][l][4][2],Y[k][l][4][3]]):
> X[k][3*l][3]:=array(1..3,[Y[k][l][5][1],Y[k][l][5][2],Y[k][l][5][3]]):
>     od:
>     od:
>     vals:={}:
>     for q from 1 to 3^(2*d) do
>     vals:=vals union {map(evalf,X[2*d+1][q][1])}:
>     vals:=vals union {map(evalf,X[2*d+1][q][2])}:
>     vals:=vals union {map(evalf,X[2*d+1][q][3])}:
>     od:
>     with(plots):
>     pointplot3d(vals,axes=BOXED);
> end: