OE:=proc(m)    local theta, i, j, a, b, c, l, X1, X2, Y1, Y2, k, vals, p:
>     theta[0]:=m:
>     for i from 1 to 11 do
>     theta[i]:=theta[i-1]*(5-theta[i-1])
>     od:
> a[3]:=1;
>     b[3]:=1-(theta[6])*(1-product( (1+4/(2-theta[n+6])), n=0..5))/4:    
>     c[3]:=b[3]:
>     for j from 2 to 0 by -1 do
>      a[j]:=((4-theta[2*j+1])*(a[j+1]+b[j+1])+2*c[j+1])/((2-theta[2*j+1])*(5-theta[2*j+1])):
> l[j]:=((4-theta[2*j+1])*(c[j+1]+b[j+1])+2*a[j+1])/((2-theta[2*j+1])*(5-theta[2*j+1])):
> b[j]:=((4-theta[2*j])*(a[j]+b[j+1])+2*l[j])/((2-theta[2*j])*(5-theta[2*j])):
> c[j]:=((4-theta[2*j])*(a[j]+l[j])+2*b[j+1])/((2-theta[2*j])*(5-theta[2*j])):
>     od:
>     for i from -1 to -5 by -1 do
>     theta[i]:=(5-sqrt(25-4*theta[i+1]))/2:
>     od:
>     X1:=array(1..6): X2:=array(1..6): Y1:=array(2..6): Y2:=array(2..6):
>     X1[1]:=array(1..1): X1[1][1]:=array(1..3): X1[1][1][1]:=array(1..3,[0.,0.,1.]): X1[1][1][2]:=array(1..3,[1.,0.,b[0]/a[0]]*1.): X1[1][1][3]:=array(1..3,[1/2.,sqrt(3)/2.,c[0]/a[0]*1.]):
>     X2[1]:=array(1..1): X2[1][1]:=array(1..3): X2[1][1][1]:=b[0]/a[0]*1.: X2[1][1][2]:=b[1]/a[0]*1.: X2[1][1][3]:=((4.-theta[0])*(l[0]+b[1])+2*a[0])/((2.-theta[0])*(5.-theta[0])):
>     for k from 2 to 6 do
>     X1[k]:=array(1..(3^(k-1))): Y1[k]:=array(1..(3^(k-2))):X2[k]:=array(1..(3^(k-1))): Y2[k]:=array(1..(3^(k-2))):
>     for l from 1 to (3^(k-2)) do
>     X1[k][3*l-2]:=array(1..3): X1[k][3*l-1]:=array(1..3): X1[k][3*l]:=array(1..3): Y1[k][l]:=array(1..6): 
>     Y1[k][l][1]:=array(1..3,[X1[k-1][l][1][1],X1[k-1][l][1][2],X1[k-1][l][1][3]]): Y1[k][l][3]:=array(1..3,[X1[k-1][l][2][1],X1[k-1][l][2][2],X1[k-1][l][2][3]]): Y1[k][l][5]:=array(1..3,[X1[k-1][l][3][1],X1[k-1][l][3][2],X1[k-1][l][3][3]]): Y1[k][l][2]:=array(1..3,[(X1[k-1][l][1][1]+X1[k-1][l][2][1])/2,X1[k-1][l][1][2],((4.-theta[1-k])*(X1[k-1][l][1][3]+X1[k-1][l][2][3])+2*X1[k-1][l][3][3])/((2.-theta[1-k])*(5.-theta[1-k]))]): Y1[k][l][4]:=array(1..3,[(X1[k-1][l][2][1]+X1[k-1][l][3][1])/2,(X1[k-1][l][2][2]+X1[k-1][l][3][2])/2,((4.-theta[1-k])*(X1[k-1][l][2][3]+X1[k-1][l][3][3])+2*X1[k-1][l][1][3])/((2.-theta[1-k])*(5.-theta[1-k]))]): Y1[k][l][6]:=array(1..3,[(X1[k-1][l][3][1]+X1[k-1][l][1][1])/2,(X1[k-1][l][3][2]+X1[k-1][l][1][2])/2,((4-theta[1-k])*(X1[k-1][l][3][3]+X1[k-1][l][1][3])+2*X1[k-1][l][2][3])/((2.-theta[1-k])*(5.-theta[1-k]))]):
>     X1[k][3*l-2]:=array(1..3): X1[k][3*l-2][1]:=array(1..3,[Y1[k][l][1][1],Y1[k][l][1][2],Y1[k][l][1][3]]): X1[k][3*l-2][2]:=array(1..3,[Y1[k][l][2][1],Y1[k][l][2][2],Y1[k][l][2][3]]): X1[k][3*l-2][3]:=array(1..3,[Y1[k][l][6][1],Y1[k][l][6][2],Y1[k][l][6][3]]):     X1[k][3*l-1]:=array(1..3): X1[k][3*l-1][1]:=array(1..3,[Y1[k][l][2][1],Y1[k][l][2][2],Y1[k][l][2][3]]): X1[k][3*l-1][2]:=array(1..3,[Y1[k][l][3][1],Y1[k][l][3][2],Y1[k][l][3][3]]): X1[k][3*l-1][3]:=array(1..3,[Y1[k][l][4][1],Y1[k][l][4][2],Y1[k][l][4][3]]): X1[k][3*l]:=array(1..3): X1[k][3*l][1]:=array(1..3,[Y1[k][l][6][1],Y1[k][l][6][2],Y1[k][l][6][3]]): X1[k][3*l][2]:=array(1..3,[Y1[k][l][4][1],Y1[k][l][4][2],Y1[k][l][4][3]]):
> X1[k][3*l][3]:=array(1..3,[Y1[k][l][5][1],Y1[k][l][5][2],Y1[k][l][5][3]]): 
>     Y2[k][l][1]:=X2[k-1][l][1]: Y2[k][l][3]:=X2[k-1][l][2]: Y2[k][l][5]:=X2[k-1][l][3]: Y2[k][l][2]:=((4.-theta[1-k])*(X2[k-1][l][1]+X2[k-1][l][2])+2*X2[k-1][l][3])/((2.-theta[1-k])*(5.-theta[1-k])): Y2[k][l][4]:=((4.-theta[1-k])*(X2[k-1][l][2]+X2[k-1][l][3])+2*X2[k-1][l][1])/((2.-theta[1-k])*(5.-theta[1-k])): Y2[k][l][6]:=((4.-theta[1-k])*(X2[k-1][l][3]+X2[k-1][l][1])+2*X2[k-1][l][2])/((2.-theta[1-k])*(5.-theta[1-k])):
>     X2[k][3*l-2]:=array(1..3): X2[k][3*l-2][1]:=Y2[k][l][1]: X2[k][3*l-2][2]:=Y2[k][l][2]: X2[k][3*l-2][3]:=Y2[k][l][6]:     X2[k][3*l-1]:=array(1..3): X2[k][3*l-1][1]:=Y2[k][l][2]: X2[k][3*l-1][2]:=Y2[k][l][3]: X2[k][3*l-1][3]:=Y2[k][l][4]: X2[k][3*l]:=array(1..3): X2[k][3*l][1]:=Y2[k][l][6]: X2[k][3*l][2]:=Y2[k][l][4]:
> X2[k][3*l][3]:=Y2[k][l][5]:
>     
>     od:
>     od:
>     vals:={}:
>     for p from 1 to 9 do
>     vals:=vals union {array(1..3,[X1[3][p][1][1],X1[3][p][1][2],X2[3][p][1]+X1[3][p][1][3]/theta[0]])}:
>     vals:=vals union {array(1..3,[X1[3][p][2][1],X1[3][p][2][2],X2[3][p][2]+X1[3][p][2][3]/theta[0]])}:
>     vals:=vals union {array(1..3,[X1[3][p][3][1],X1[3][p][3][2],X2[3][p][3]+X1[3][p][3][3]/theta[0]])}:
>     od:
>     vals:=map(evalf,vals):
>     print(vals);
>     vals:={}:
>     for p from 1 to 243 do
>     vals:=vals union {array(1..3,[X1[6][p][1][1],X1[6][p][1][2],X2[6][p][1]+X1[6][p][1][3]/theta[0]])}:
>     vals:=vals union {array(1..3,[X1[6][p][2][1],X1[6][p][2][2],X2[6][p][2]+X1[6][p][2][3]/theta[0]])}:
>     vals:=vals union {array(1..3,[X1[6][p][3][1],X1[6][p][3][2],X2[6][p][3]+X1[6][p][3][3]/theta[0]])}:
>     od:
>     vals:=map(evalf,vals):
>     with(plots):
>     pointplot3d(vals,axes=BOXED,title="Lambda0="||m);
> end: