Useful functions/definitions:

(* Numerically check if a square matrix is the identity matrix to within some accuracy *)
isIdentityMatrix[r_]:=(Flatten[N[r-IdentityMatrix[Length[r]]]].Flatten[N[r-IdentityMatrix[Length[r]]]]<0.001);

(* Numerically check if any matrix is the zero matrix to within some accuracy *)
isZeroMatrix[r_]:=(Flatten[N[r]].Flatten[N[r]]<0.001);

(* Numerically test and remove duplicate matrices from a list *)
reduce[list_]:=Union[list,SameTest->(isZeroMatrix[#1-#2]&)];

(* Given a list of square matrices g, this function combines all possible pairs of matrices in all possible orders, and numerically removes all duplicate matrices *)
iterate[g_]:=reduce[FullSimplify[#1.#2]&@@@Tuples[g,2]];

(* Given a generating set of matrices, this function creates a group. If the generators aren't actually a generating set, 
this code runs forever! *)
makeGroup[generators_]:=FixedPoint[iterate,Join[generators,{IdentityMatrix[3]}]];

(* Find what integer power you have to raise a matrix to to get the identity. eg,
given integer "m", the command getMatrixOrder[RotationMatrix[2Pi/m,{1,0,0}]] will return "m". *)
getMatrixOrder[r_]:=Module[{n,mat},n=1;mat=r;
While[n<100&& Not[isIdentity[mat]],
mat=N[mat.r];
n=n+1;
];
n
];

(* Find the invariant axis of a three dimensional rotation matrix by selecting for the  only eigenvector with eigenvalue one. *)
rotAxis[r_]:=Select[Transpose[Eigensystem[r]],#[[1]]==1&][[1,2]];

(* In finding axes of rotation, we wish to remove all duplicate axes.
An axis is a line passing through the origin, and I choose to represent
it by a unit vector. However, the unit vector {-1,0,0} and {1,0,0} both 
represent the same axis -- the line {x,0,0}. We can negate the unit 
vector arbitrarily, and we can use this freedom to ensure that the leftmost number 
is positive. This function does just that: negates a 3D vector to ensure the 
leftmost component is positive. *)
orderAxis[v_]:=If[v[[1]]<0,-v,If[v[[1]]==0,
If[v[[2]]<0,-v,If[v[[2]]==0,
If[v[[3]]<0,-v,v],v]],v]];

normalize[v_]:=Simplify[v/Sqrt[v.v]];

getRotationAxes[group_]:=Module[{lists,myfunction},
(* "lists" is of the form {{order,axis},{order,axis},{order,axis}}
where "axis" is a vector whose leftmost coordinate is positive, 
and "order" is an integer between 2 and 5 inclusive. *)
lists={getMatrixOrder[#],orderAxis[rotAxis[#]]}&/@group;
(* Select elements from list with a given order, remove duplicates, and normalize *)
myfunction[n_]:=normalize/@reduce[Select[lists,#[[1]]==n&][[All,2]]];
{{2,myfunction[2]} (* all axes of order two *),
{3,myfunction[3]} (* all axes of order three *),
{4,myfunction[4]} (* all axes of order four *),
{5,myfunction[5]}}];

orbit[point_,group_]:=reduce[#.point&/@group];

The Groups and their Axes of Rotation

(* Tetrahedron generators, group elements, and all axes of rotation *)
tetGenerators={RotationMatrix[2 Pi/3,{1,0,-1/Sqrt[2]}],RotationMatrix[2 Pi/3,{0,1,1/Sqrt[2]}]};
tetG=makeGroup[tetGenerators];
tetaxes=getRotationAxes[tetG];
tetaxes2=tetaxes[[1,2]];
tetaxes3=tetaxes[[2,2]];

(* Octahedron generators, group elements, and all axes of rotation *)
octGenerators={RotationMatrix[Pi/2,{1,0,0}],RotationMatrix[Pi/2,{0,1,0}]};
octG=makeGroup[octGenerators];
octaxes=getRotationAxes[octG];
octaxes2=octaxes[[1,2]];
octaxes3=octaxes[[2,2]];
octaxes4=octaxes[[3,2]];

(* Icosahedron generators, group elements, and all axes of rotation *)
f=(1+Sqrt[5])/2;
icoGenerators={{{-1,0,0},{0,-1,0},{0,0,01}},{{0,0,1},{1,0,0},{0,1,0}},1/2 {{1,-f,1/f},{f,1/f,-1},{1/f,1,f}}};
icoG=makeGroup[icoGenerators];
icoaxes=getRotationAxes[icoG];
icoaxes2=icoaxes[[1,2]];
icoaxes3=icoaxes[[2,2]];
icoaxes5=icoaxes[[4,2]];

tetV=orbit[tetaxes3[[2]],tetG]; (* tetrahedron vertices *)
cubeoctV=orbit[octaxes2[[2]],octG];(* cuboctahedron vertices *)
cubeV=orbit[octaxes3[[1]],octG]; (* cube vertices *)
octV=orbit[octaxes2[[1]],octG]; (* octahedron vertices *)
icosadodecaV=orbit[icoaxes2[[1]],icoG]; (* icosadodecahedron vertices *)
dodecaV=orbit[icoaxes3[[1]],icoG]; (* dodecahedron vertices *)
icosaV=orbit[icoaxes5[[1]],icoG]; (* icosahedron vertices *)

Some nice drawings:

drawhull[verts_]:=Graphics3D[{EdgeForm[],GraphicsComplex[MeshCoordinates@#,MeshCells[#,2]]&[ConvexHullMesh[verts]]},Boxed->False];
drawhull/@{tetV,cubeoctV,cubeV,octV,icosadodecaV,dodecaV,icosaV}

drawhullAxes[verts_,axes_]:=Graphics3D[{Thick,Transpose[{{Blue,Green,Yellow,Red},(Line[{1.2#,-1.2#}]&/@#)&/@axes[[All,2]]}],EdgeForm[],GraphicsComplex[MeshCoordinates@#,MeshCells[#,2]]&[ConvexHullMesh[verts]]},Boxed->False];
drawhullAxes@@@{{tetV,tetaxes},{cubeoctV,octaxes},{cubeV,octaxes},{octV,octaxes},{icosadodecaV,icoaxes},{dodecaV,icoaxes},{icosaV,icoaxes}}

Text Output:

jsstring[list_]:=StringReplace[ToString[Round[Chop[N[list]],.00001]],{"{"->"[","}"->"]","*^"->"e"}];
"var tetMatrices="<>jsstring[tetG]<>";\n\n"<>
"var tetAxes2="<>jsstring[tetaxes2]<>";\n\n"<>
"var tetAxes3="<>jsstring[tetaxes3]<>";\n\n"<>
"var octMatrices="<>jsstring[octG]<>";\n\n"<>
"var octAxes2="<>jsstring[octaxes2]<>";\n\n"<>
"var octAxes3="<>jsstring[octaxes3]<>";\n\n"<>
"var octAxes4="<>jsstring[octaxes4]<>";\n\n"<>
"var icoMatrices="<>jsstring[icoG]<>";\n\n"<>
"var icoAxes="<>jsstring[icoaxes2]<>";\n\n"<>
"var icoAxes="<>jsstring[icoaxes3]<>";\n\n"<>
"var icoAxes="<>jsstring[icoaxes5]<>";\n\n"<>
"var tetV="<>jsstring[tetV]<>";\n\n"<>
"var cubeoctV="<>jsstring[cubeoctV]<>";\n\n"<>
"var cubeV="<>jsstring[cubeV]<>";\n\n"<>
"var octV="<>jsstring[octV]<>";\n\n"<>
"var icosadodecaV="<>jsstring[icosadodecaV]<>";\n\n"<>
"var dodecaV="<>jsstring[dodecaV]<>";\n\n"<>
"var icosaV="<>jsstring[icosaV]<>";\n\n"