Consider a cube with a sphere inside it:
We can fit a smaller sphere snugly into each of the corners; finding the radius involves solving a quadratic equation, and the answer is X = Y = Z = T0 = sqrt(3)-1, radius = R0 = 2-sqrt(3).
The next stage is more interesting. We want to maximise the function min(1-X, 1-Y, 1-Z, sqrt(X^2+Y^2+Z^2)-1, sqrt((X-T0)^2+(Y-T0)^2+(Z-T0)^2)-R0^2, the radius of the largest ball that would fit at position (X,Y,Z). Unsurprisingly, this objective function has a fair number of local maxima. First, I tried using Excel's Solver function, which found me another point lying along the diagonal.
But Excel Solver gives approximate answers and doesn't seem very good at walking along gradients. Time to switch to a heftier tool: a short Python script which walks from a point X in the direction of grad(f)(X), in exponentially-increasing steps until f(X) stops increasing. Here's the script, including both a grad f method and a stupid process that samples at points in a cube and picks the best value.
The stupider method works best, as is depressingly often the case. The problem is that I'm not sure of a good way to compute grad f numerically for a function, particularly near a local minimum where it's in any case close to zero. And the function is discontinuous, so grad f is not well-defined (which I realised by noticing that my grad function was often returning non-zero grads at points that the stupid process could not refine).
Maybe the best way is analytic. The lazy-sampling method suggests that the correct point (where "correct" means that the radius is largest) lies equidistant from a side plate, the central sphere and a corner sphere - there are 24 such points, since there are 8 corner spheres and 3 possible side plates for each sphere. But having three quantities forced equal gives only two quadratic equations in three unknowns, and having covered three pieces of paper with working to solve these, I gave up.
We know the radius to fourteen decimal places. So, using a lattice-reduction technique for finding minimal polynomials, we can observe that the radius satisfies 25r^4 - 90r^3 + 94r^2 - 36r + 4 = 0. This is, to put it mildly, unenlightening. We can check this is sane by noticing that, if the co-ordinate is [1-r, 1-r, v], we have v = sqrt(6r - r^2 - 1), and that v is then defined over the same field; also, the radius-polynomial has Galois group C2 x C2 (is it not wonderful to have magma at one's fingertips), which is what you'd expect from a construction which, hopefully, could be done with ruler and compasses. If anyone can explain to me why it is natural to note that (5r^2 - 9r + 2)^2 = 7r^2, or indeed any other pair of quadratics, I'll give them a Mars bar.
OK, this is now quite pretty, but not very mathematical: the correct answer must be to figure out some elegant way of writing the position and radius of the sphere tangent to three arbitrarily-located spheres of arbitrary radius (modelling the sides of the cube as parts of the sides of six spheres of infinite radius). And then you'd just iterate down, seeing which of the triplets of known spheres give the largest radius for the tangent sphere.
Renderings for this page were done with Persistence of Vision, using this input file (it's fairly wonderful that ray-tracing is nowadays quick enough for almost interactive use — that scene renders in 19 seconds), and converted to .PNG format with GIMP; sadly the use of Magma means this is not a purely Free Software project.