**Edit:** a preprint concerning this problem can now be found on the arXiv: http://arxiv.org/abs/1407.0683

Let me give an exhaustive answer. Croft closes his paper with a list of the unsolved cases:

Here $\kappa$ denotes the maximal edge length of the inner polyhedron. The ratio of the volumes can be easily computed when $\kappa$ is known.

Let's fill in the blank spaces that "denote ignorance"!

Let's say we want to find the largest polyhedron similar to a regular polyhedron $P$ that fits into a regular polyhedron $Q$ with edge length equal to $1$. We set up the following optimization problem.

*real variables:* $e$, and for every vertex $v$ of $P$: $(x_v,y_v,z_v)$

*objective:* maximize $e$

*linear constraints:* For each vertex $v$ in $P$ and each hyperplane $h^+$ defining $Q$: $$(x_v,y_v,z_v)\in h^+.$$

*quadratic constraints:* For each edge $(v,w)$ in $P$: $$e=(x_v-x_w)^2+(y_v-y_w)^2+(z_v-z_w)^2.$$
If $Q=C$ or $Q=D$ we have to include more constraints to make sure that the coordinates $(x_v,y_v,z_v)$ really give a polyhedron similar to $P$.

In other words: find coordinates of a polyhedron similar to $P$, which are contained in $Q$ with the square of the edge length being equal to $e$. Maximize $e$ to find the solution.

This non-linear programming problem with quadratic constraints can be effectively solved with SCIP, which uses branch-and-bound methods. (Many thanks to Ambros Gleixner for help with SCIP!). In order to run fast enough one has to set up the problem a bit more clever and remove redundant constraints.
I ran the program to solve the unsolved cases and also checked to known ones. Here is a table with the results:

I emphasized new results (where Croft had blanks). Notice also that my values differ from Croft's in three places, which I marked with an exclamation mark.

- In the case $T$ in $D$ the number is off by a power of ten.
- In the case $O$ in $D$ the rounded decimal expansion should read $1.851$, not $1.1850$, since the number really is $1.85122958682191611960\ldots$.
- In the dual case $D$ in $C$ the correct decimal expansion reads $0.3942834797251374518168\ldots,$ so the rounded number should be $0.394$ not $0.395$.

However Croft gives correct exact results ($\tau$ denotes the golden ratio) and only makes these minor mistakes in the decimal approximation.

Using similar arguments about immobility as Croft, we can calculate exact results for the open pairs. Let me give you exact values for two open pairs.

For $D$ in $I$ the maximum is $$\frac{15-\sqrt{5}}{22}\approx 0.58017872829546410470867392\ldots.$$
For $C$ in $I$ the maximum is $$\frac{5+7\sqrt{5}}{22}\approx 0.93874890193175126703928253\ldots.$$

In some cases it is harder to find the algebraic values, but in all cases it can be done, see the arXiv preprint above for details.

I used sage to generate a few images. Click here to view a 3d animation.

The same optimization ansatz will also work in other dimensions and even for non-regular polyhedra.