A3:=Matrix([[1,1,1,6,6,6,3,3,3],
 [1,1,1,6,6,6,3,3,3],
 [1,1,1,6,6,6,3,3,3],
 [2,2,2,5,5,5,4,4,4],
 [2,2,2,5,5,5,4,4,4],
 [2,2,2,5,5,5,4,4,4]]):


rubik3_squares := proc () 
local i; 
for i from 1 to 4 do 
F1||i:=[3,-3+2*(i-1),3];
F2||i:=[3,-3+2*(i-1),1];
F3||i:=[3,-3+2*(i-1),-1];
F4||i:=[3,-3+2*(i-1),-3];
R1||i:=[-3+2*(i-1),3,3];
R2||i:=[-3+2*(i-1),3,1];
R3||i:=[-3+2*(i-1),3,-1];
R4||i:=[-3+2*(i-1),3,-3];
B1||i:=[-3,-3+2*(i-1),3];
B2||i:=[-3,-3+2*(i-1),1];
B3||i:=[-3,-3+2*(i-1),-1];
B4||i:=[-3,-3+2*(i-1),-3];
L1||i:=[-3+2*(i-1),-3,3];
L2||i:=[-3+2*(i-1),-3,1];
L3||i:=[-3+2*(i-1),-3,-1];
L4||i:=[-3+2*(i-1),-3,-3];
U1||i:=[-3+2*(i-1),3,3];
U2||i:=[-3+2*(i-1),1,3];
U3||i:=[-3+2*(i-1),-1,3];
U4||i:=[-3+2*(i-1),-3,3];
D1||i:=[-3+2*(i-1),3,-3];
D2||i:=[-3+2*(i-1),1,-3];
D3||i:=[-3+2*(i-1),-1,-3];
D4||i:=[-3+2*(i-1),-3,-3];
end do;
end proc:

cube3_draw := proc (A) 
local i, j, help, 
face_F11,face_F12,face_F13,
face_F21,face_F22,face_F23,
face_F31,face_F32,face_F33,
face_B11,face_B12,face_B13,
face_B21,face_B22,face_B23,
face_B31,face_B32,face_B33,
face_U11,face_U12,face_U13,
face_U21,face_U22,face_U23,
face_U31,face_U32,face_U33,
face_R11,face_R12,face_R13,
face_R21,face_R22,face_R23,
face_R31,face_R32,face_R33,
face_L11,face_L12,face_L13,
face_L21,face_L22,face_L23,
face_L31,face_L32,face_L33, 
face_D11,face_D12,face_D13,
face_D21,face_D22,face_D23,
face_D31,face_D32,face_D33, L; 
rubik3_squares();
help:=[red,blue,coral,plum,green,yellow];
for i from 1 to 6 do
 color||i:=help[i];
end do;
face_F11:=plots[polygonplot3d]([F11,F12,F22,F21],style=patch,
color=color||(A[1,1]));
face_F12:=plots[polygonplot3d]([F12,F13,F23,F22],style=patch,
color=color||(A[1,2]));
face_F13:=plots[polygonplot3d]([F13,F14,F24,F23],style=patch,
color=color||(A[1,3]));
face_F21:=plots[polygonplot3d]([F21,F22,F32,F31],style=patch,
color=color||(A[2,1]));
face_F22:=plots[polygonplot3d]([F22,F23,F33,F32],style=patch,
color=color||(A[2,2]));
face_F23:=plots[polygonplot3d]([F23,F24,F34,F33],style=patch,
color=color||(A[2,3]));
face_F31:=plots[polygonplot3d]([F31,F32,F42,F41],style=patch,
color=color||(A[3,1]));
face_F32:=plots[polygonplot3d]([F32,F33,F43,F42],style=patch,
color=color||(A[3,2]));
face_F33:=plots[polygonplot3d]([F33,F34,F44,F43],style=patch,
color=color||(A[3,3]));
face_B11:=plots[polygonplot3d]([B11,B12,B22,B21],style=patch,
color=color||(A[1,7]));
face_B12:=plots[polygonplot3d]([B12,B13,B23,B22],style=patch,
color=color||(A[1,8]));
face_B13:=plots[polygonplot3d]([B13,B14,B24,B23],style=patch,
color=color||(A[1,9]));
face_B21:=plots[polygonplot3d]([B21,B22,B32,B31],style=patch,
color=color||(A[2,7]));
face_B22:=plots[polygonplot3d]([B22,B23,B33,B32],style=patch,
color=color||(A[2,8]));
face_B23:=plots[polygonplot3d]([B23,B24,B34,B33],style=patch,
color=color||(A[2,9]));
face_B31:=plots[polygonplot3d]([B31,B32,B42,B41],style=patch,
color=color||(A[3,7]));
face_B32:=plots[polygonplot3d]([B32,B33,B43,B42],style=patch,
color=color||(A[3,8]));
face_B33:=plots[polygonplot3d]([B33,B34,B44,B43],style=patch,
color=color||(A[3,9]));
face_U11:=plots[polygonplot3d]([U11,U12,U22,U21],style=patch,
color=color||(A[4,1]));
face_U21:=plots[polygonplot3d]([U12,U13,U23,U22],style=patch,
color=color||(A[4,2]));
face_U31:=plots[polygonplot3d]([U13,U14,U24,U23],style=patch,
color=color||(A[4,3]));
face_U12:=plots[polygonplot3d]([U21,U22,U32,U31],style=patch,
color=color||(A[5,1]));
face_U22:=plots[polygonplot3d]([U22,U23,U33,U32],style=patch,
color=color||(A[5,2]));
face_U32:=plots[polygonplot3d]([U23,U24,U34,U33],style=patch,
color=color||(A[5,3]));
face_U13:=plots[polygonplot3d]([U31,U32,U42,U41],style=patch,
color=color||(A[6,1]));
face_U23:=plots[polygonplot3d]([U32,U33,U43,U42],style=patch,
color=color||(A[6,2]));
face_U33:=plots[polygonplot3d]([U33,U34,U44,U43],style=patch,
color=color||(A[6,3]));
face_R11:=plots[polygonplot3d]([R11,R12,R22,R21],style=patch,
color=color||(A[1,4]));
face_R12:=plots[polygonplot3d]([R12,R13,R23,R22],style=patch,
color=color||(A[1,5]));
face_R13:=plots[polygonplot3d]([R13,R14,R24,R23],style=patch,
color=color||(A[1,6]));
face_R21:=plots[polygonplot3d]([R21,R22,R32,R31],style=patch,
color=color||(A[2,4]));
face_R22:=plots[polygonplot3d]([R22,R23,R33,R32],style=patch,
color=color||(A[2,5]));
face_R23:=plots[polygonplot3d]([R23,R24,R34,R33],style=patch,
color=color||(A[2,6]));
face_R31:=plots[polygonplot3d]([R31,R32,R42,R41],style=patch,
color=color||(A[3,4]));
face_R32:=plots[polygonplot3d]([R32,R33,R43,R42],style=patch,
color=color||(A[3,5]));
face_R33:=plots[polygonplot3d]([R33,R34,R44,R43],style=patch,
color=color||(A[3,6]));
face_L11:=plots[polygonplot3d]([L11,L12,L22,L21],style=patch,
color=color||(A[4,4]));
face_L12:=plots[polygonplot3d]([L12,L13,L23,L22],style=patch,
color=color||(A[4,5]));
face_L13:=plots[polygonplot3d]([L13,L14,L24,L23],style=patch,
color=color||(A[4,6]));
face_L21:=plots[polygonplot3d]([L21,L22,L32,L31],style=patch,
color=color||(A[5,4]));
face_L22:=plots[polygonplot3d]([L22,L23,L33,L32],style=patch,
color=color||(A[5,5]));
face_L23:=plots[polygonplot3d]([L23,L24,L34,L33],style=patch,
color=color||(A[5,6]));
face_L31:=plots[polygonplot3d]([L31,L32,L42,L41],style=patch,
color=color||(A[6,4]));
face_L32:=plots[polygonplot3d]([L32,L33,L43,L42],style=patch,
color=color||(A[6,5]));
face_L33:=plots[polygonplot3d]([L33,L34,L44,L43],style=patch,
color=color||(A[6,6]));
face_D11:=plots[polygonplot3d]([D11,D12,D22,D21],style=patch,
color=color||(A[4,7]));
face_D21:=plots[polygonplot3d]([D12,D13,D23,D22],style=patch,
color=color||(A[4,8]));
face_D31:=plots[polygonplot3d]([D13,D14,D24,D23],style=patch,
color=color||(A[4,9]));
face_D12:=plots[polygonplot3d]([D21,D22,D32,D31],style=patch,
color=color||(A[5,7]));
face_D22:=plots[polygonplot3d]([D22,D23,D33,D32],style=patch,
color=color||(A[5,8]));
face_D32:=plots[polygonplot3d]([D23,D24,D34,D33],style=patch,
color=color||(A[5,9]));
face_D13:=plots[polygonplot3d]([D31,D32,D42,D41],style=patch,
color=color||(A[6,7]));
face_D23:=plots[polygonplot3d]([D32,D33,D43,D42],style=patch,
color=color||(A[6,8]));
face_D33:=plots[polygonplot3d]([D33,D34,D44,D43],style=patch,
color=color||(A[6,9]));
L:= [face_F11,face_F12,face_F13,
face_F21,face_F22,face_F23,
face_F31,face_F32,face_F33,
face_B11,face_B12,face_B13,
face_B21,face_B22,face_B23,
face_B31,face_B32,face_B33,
face_U11,face_U12,face_U13,
face_U21,face_U22,face_U23,
face_U31,face_U32,face_U33,
face_R11,face_R12,face_R13,
face_R21,face_R22,face_R23,
face_R31,face_R32,face_R33,
face_L11,face_L12,face_L13,
face_L21,face_L22,face_L23,
face_L31,face_L32,face_L33,
face_D11,face_D12,face_D13,
face_D21,face_D22,face_D23,
face_D31,face_D32,face_D33]; 
return L 
end proc:

move3_r := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[1,3],[2,3],[3,3],[4,3],[4,2],[4,1],[1,9],
[2,9],[3,9],[4,9],[4,8],[4,7],[1,4],[1,6],[3,6],
[3,4],[1,5],[2,6],[3,5],[2,4]};
AA[4,3]:=A[3,3]; 
AA[4,2]:=A[2,3];
AA[4,1]:=A[1,3]; #front to up
AA[1,9]:=A[4,3];
AA[2,9]:=A[4,2]; 
AA[3,9]:=A[4,1]; #up to back
AA[4,9]:=A[3,9];
AA[4,8]:=A[2,9]; 
AA[4,7]:=A[1,9]; #back to down
AA[3,3]:=A[4,7];
AA[2,3]:=A[4,8]; 
AA[1,3]:=A[4,9]; #up to back
AA[1,4]:=A[1,6];
AA[1,6]:=A[3,6];
AA[3,6]:=A[3,4];
AA[3,4]:=A[1,4];
AA[1,5]:=A[2,6];
AA[2,6]:=A[3,5];
AA[3,5]:=A[2,4];
AA[2,4]:=A[1,5];
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_r

move3_u := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[1,1],[1,2],[1,3],[1,4],[1,5],[1,6],[1,7],
[1,8],[1,9],[4,4],[4,5],[4,6],[6,1],[6,3],[4,3],
[4,1],[6,2],[5,3],[4,2],[5,1]};
AA[1,6]:=A[1,9]; 
AA[1,5]:=A[1,8];
AA[1,4]:=A[1,7]; #front to right
AA[1,9]:=A[4,4];
AA[1,8]:=A[4,5]; 
AA[1,7]:=A[4,6]; #back to right
AA[4,6]:=A[1,3];
AA[4,5]:=A[1,2]; 
AA[4,4]:=A[1,1]; #back to down
AA[1,1]:=A[1,6];
AA[1,2]:=A[1,5]; 
AA[1,3]:=A[1,4]; #up to back
AA[4,1]:=A[6,1];
AA[6,1]:=A[6,3];
AA[6,3]:=A[4,3];
AA[4,3]:=A[4,1];
AA[6,2]:=A[5,3];
AA[5,3]:=A[4,2];
AA[4,2]:=A[5,1];
AA[5,1]:=A[6,2];
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_u

move3_l := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[1,1],[2,1],[3,1],[1,7],[2,7],[3,7],[6,7],
[6,8],[6,9],[6,1],[6,2],[6,3],[4,4],[6,4],[6,6],[4,6],
[4,5],[5,4],[6,5],[5,6]};
AA[1,1]:=A[6,1]; 
AA[2,1]:=A[6,2];
AA[3,1]:=A[6,3]; #front to up
AA[6,7]:=A[3,1];
AA[6,8]:=A[2,1]; 
AA[6,9]:=A[1,1]; #back to up
AA[1,7]:=A[6,7];
AA[2,7]:=A[6,8]; 
AA[3,7]:=A[6,9]; #back to down
AA[6,1]:=A[3,7];
AA[6,2]:=A[2,7]; 
AA[6,3]:=A[1,7]; #up to back
AA[6,4]:=A[6,6];
AA[4,4]:=A[6,4];
AA[4,6]:=A[4,4];
AA[6,6]:=A[4,6];
AA[5,4]:=A[6,5];
AA[6,5]:=A[5,6];
AA[5,6]:=A[4,5];
AA[4,5]:=A[5,4];
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_l

move3_f := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[1,6],[2,6],[3,6],[4,9],[5,9],[6,9],[6,6],
[5,6],[4,6],[6,3],[5,3],[4,3],[1,1],[3,1],
[3,3],[1,3],[1,2],[2,1],[3,2],[2,3]};
AA[1,1]:=A[3,1];
AA[3,1]:=A[3,3];
AA[3,3]:=A[1,3];
AA[1,3]:=A[1,1];
AA[1,2]:=A[2,1];
AA[2,1]:=A[3,2];
AA[3,2]:=A[2,3];
AA[2,3]:=A[1,2];
AA[1,6]:=A[6,3]; 
AA[2,6]:=A[5,3];
AA[3,6]:=A[4,3]; #up to rear
AA[4,9]:=A[1,6];
AA[5,9]:=A[2,6]; 
AA[6,9]:=A[3,6]; #rear to down
AA[6,6]:=A[4,9];
AA[5,6]:=A[5,9]; 
AA[4,6]:=A[6,9]; #down to left
AA[6,3]:=A[6,6];
AA[5,3]:=A[5,6]; 
AA[4,3]:=A[4,6]; #left to up
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_f

move3_b := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[4,1],[5,1],[6,1],[6,4],[5,4],
[4,4],[4,7],[5,7],[6,7],[3,4],
[2,4],[1,4],[1,8],[1,9],[1,7],
[2,9],[2,7],[3,7],[3,9],[3,8]};
AA[1,8]:=A[2,9];
AA[2,9]:=A[3,8];
AA[3,8]:=A[2,7];
AA[2,7]:=A[1,8];
AA[1,9]:=A[3,9];
AA[3,9]:=A[3,7];
AA[3,7]:=A[1,7];
AA[1,7]:=A[1,9];
AA[4,1]:=A[3,4]; 
AA[5,1]:=A[2,4];
AA[6,1]:=A[1,4]; #right to up
AA[6,4]:=A[6,1];
AA[5,4]:=A[5,1]; 
AA[4,4]:=A[4,1]; #up to left
AA[4,7]:=A[6,4];
AA[5,7]:=A[5,4]; 
AA[6,7]:=A[4,4]; #down to left
AA[3,4]:=A[6,7];
AA[2,4]:=A[5,7]; 
AA[1,4]:=A[4,7]; #left to up
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_b

move3_d := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[3,6],[3,5],[3,4],[3,9],[3,8],[3,7],[6,4],
[6,5],[6,6],[3,1],[3,2],[3,3],[6,9],[6,7],[4,7],[4,9],
[5,9],[6,8],[5,7],[4,8]};
AA[6,9]:=A[6,7];
AA[6,7]:=A[4,7];
AA[4,7]:=A[4,9];
AA[4,9]:=A[6,9];
AA[5,9]:=A[6,8];
AA[6,8]:=A[5,7];
AA[5,7]:=A[4,8];
AA[4,8]:=A[5,9];
AA[3,6]:=A[3,1]; 
AA[3,5]:=A[3,2];
AA[3,4]:=A[3,3]; #front to right
AA[3,9]:=A[3,6];
AA[3,8]:=A[3,5]; 
AA[3,7]:=A[3,4]; #up to left
AA[6,4]:=A[3,9];
AA[6,5]:=A[3,8]; 
AA[6,6]:=A[3,7]; #down to left
AA[3,1]:=A[6,4];
AA[3,2]:=A[6,5]; 
AA[3,3]:=A[6,6]; #left to up
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_d

move3_mr := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[1,2],[2,2],[3,2],[5,7],[5,8],[5,9],[1,8],[2,8],[3,8],[5,1],[5,2],[5,3]};
AA[1,2]:=A[5,9];
AA[5,1]:=A[1,2];
AA[3,8]:=A[5,1];
AA[5,9]:=A[3,8];
AA[2,2]:=A[5,8];
AA[5,2]:=A[2,2];
AA[2,8]:=A[5,2];
AA[5,8]:=A[2,8];
AA[3,2]:=A[5,7]; 
AA[5,3]:=A[3,2];
AA[1,8]:=A[5,3];
AA[5,7]:=A[1,8];
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_d

move3_md := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[2,1],[2,2],[2,3],[2,4],[2,5],[2,6],[2,9],[2,8],[2,7],[5,4],[5,5],[5,6]};
AA[2,2]:=A[5,5];
AA[2,5]:=A[2,2];
AA[2,8]:=A[2,5];
AA[5,5]:=A[2,8];
AA[2,1]:=A[5,4];
AA[2,6]:=A[2,1];
AA[2,9]:=A[2,6];
AA[5,4]:=A[2,9];
AA[2,3]:=A[5,6]; 
AA[2,4]:=A[2,3];
AA[2,7]:=A[2,4];
AA[5,6]:=A[2,7];
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_d

move3_mf := proc (A) 
local j, k, AA,indices; 
AA := Matrix(6,9); 
indices:={[4,5],[5,5],[6,5],[6,2],[5,2],[4,2],[1,5],[2,5],[3,5],[6,8],[5,8],[4,8]};
AA[2,5]:=A[5,2];
AA[5,2]:=A[5,5];
AA[5,5]:=A[5,8];
AA[5,8]:=A[2,5];
AA[6,8]:=A[3,5];
AA[3,5]:=A[4,2];
AA[4,2]:=A[4,5];
AA[4,5]:=A[6,8];
AA[6,2]:=A[6,5]; 
AA[1,5]:=A[6,2];
AA[4,8]:=A[1,5];
AA[6,5]:=A[4,8];
for j from 1 to 6 do
 for k from 1 to 9 do
  if not(member([j,k],indices)) then
   AA[j,k]:=A[j,k] 
  end if;
 end do;
end do;
return AA
end proc: #move_d

move3_power:=proc(g,A)
local i,j,AA,generators,gen,power;
generators:={R,U,F,D,L,B,Mr,Md,Mf};
AA:=A;
if type(g,`^`) then 
  gen:=op(1,g);
  power:=op(2,g) mod 4
 else if member(g,generators) then 
   gen:=g;
   power:=1
 end if
end if;
if g=1 then return AA  end if;
for i from 1 to power do
 if gen=Mr then AA:=move3_mr(AA) end if;
 if gen=Md then AA:=move3_md(AA) end if;
 if gen=Mf then AA:=move3_mf(AA) end if;
 if gen=R then AA:=move3_r(AA) end if;
 if gen=L then AA:=move3_l(AA) end if;
 if gen=U then AA:=move3_u(AA) end if;
 if gen=D then AA:=move3_d(AA) end if;
 if gen=F then AA:=move3_f(AA) end if;
 if gen=B then AA:=move3_b(AA) end if;
end do;
return AA
end proc:

move3:=proc(w,A)
local i,AA,generators;
generators:={Mr,Md,Mf,R,U,F,D,L,B};
AA:=A;
if type(w,`^`) then 
  AA:=move3_power(w,AA)
end if;
if w=1 then return AA end if;
if member(w,generators) then
  AA:=move3_power(w,AA)
end if;
if type(w,list) then 
 for i from 1 to nops(w) do
   AA:=move3_power(op(i,w),AA)
 end do
end if;
return AA
end proc:

rc3 := proc (L, A, RET) local AA, LL; 
AA := move3(L,A); 
if nargs=3 then 
 if not type(RET,name) then error `wrong third argument` end if;
 RET:=AA 
end if;
LL := cube3_draw(AA); 
plots[display3d](LL,scaling=CONSTRAINED,projection=.95) 
end proc: