Determine height of box packed with spheres

Your question is like a sangaku problem. What is the distance between the centre of the large sphere and one of the smaller spheres? The large sphere centre is $\{2,2,h/2\}$, one of the small spheres is at $\{1,1,1\}$. Their separation equals the sum of their radii. That is,

Norm[{2, 2, h/2} - {1, 1, 1}] == 1 + 2

Solve the expression for $h$, and choose the positive root, $h=2(1+\sqrt{7})$.

With[{h = 2 (1 + Sqrt[7])},
   Graphics3D[{Orange,
      Sphere[{1, 1, 1}, 1], Sphere[{3, 1, 1}, 1], Sphere[{1, 3, 1}, 1],
      Sphere[{3, 3, 1}, 1], Sphere[{1, 1, h - 1}, 1], 
      Sphere[{3, 1, h - 1}, 1],
      Sphere[{1, 3, h - 1}, 1], Sphere[{3, 3, h - 1}, 1],
      Red, Sphere[{2, 2, h/2}, 2]
     }, Lighting -> "Neutral",Axes->True]]

KennyColnago's answer is good, but we're using Mathematica, so we can leave much more of the problem to it. Now, we can't avoid writing down the system of equations; for simplicity, I'll assume that the large sphere (with radius $R = 2$) is centered at the origin, and that the small spheres (with radius $r = 1$) are centered at ${x, y, z}$. The sides of the box lie at $\pm2$, and the top and bottom are at $\pm h$, so we have tangent small spheres for any combination of positive and negative $x$, $y$ and $z$. That means we know the following:

In[1]:= eqns = {
           (* each sphere is tangent to the sides of the box *)
           Abs[x] + r == 2,
           Abs[y] + r == 2,

           (* each sphere is tangent to the top of the box *)
           Abs[z] + r == h/2,

           (* because the large sphere and small sphere are tangent 
              to one another, their radii are collinear, so their 
              total length is the distance of the center of the small
              sphere from the origin. *)
           r + R == Norm[{x, y, z}],

           (* h is positive *)
           0 < h
        } /. {r -> 1, R -> 2};

Now, we can use Mathematica to solve those equations for $x$, $z$ and $h$.

In[2]:= sol = Simplify@Solve[eqns, {x, y, z, h}, Reals]
Out[2]= {{x -> -1, y -> -1, z -> -Sqrt[7], h -> 2 (1 + Sqrt[7])}, 
         {x -> -1, y -> -1, z -> Sqrt[7], h -> 2 (1 + Sqrt[7])}, 
         {x -> -1, y -> 1, z -> -Sqrt[7], h -> 2 (1 + Sqrt[7])}, 
         {x -> -1, y -> 1, z -> Sqrt[7], h -> 2 (1 + Sqrt[7])}, 
         {x -> 1, y -> -1, z -> -Sqrt[7], h -> 2 (1 + Sqrt[7])}, 
         {x -> 1, y -> -1, z -> Sqrt[7], h -> 2 (1 + Sqrt[7])}, 
         {x -> 1, y -> 1, z -> -Sqrt[7], h -> 2 (1 + Sqrt[7])},  
         {x -> 1, y -> 1, z -> Sqrt[7], h -> 2 (1 + Sqrt[7])}}

In[3]:= Length@sol
Out[3]= 8

Gratifyingly, we have eight solutions, one for each small sphere.

Now we can draw our spheres, which is easy because Mathematica has helpfully collected all the solutions as a list of rules.

In[3]:= centers = {x, y, z} /. sol
Out[3]= {{-1, -1, -Sqrt[7]}, {-1, -1, Sqrt[7]}, 
         {-1, 1, -Sqrt[7]}, {-1, 1, Sqrt[7]}, 
         {1, -1, -Sqrt[7]}, {1, -1, Sqrt[7]}, 
         {1, 1, -Sqrt[7]}, {1, 1, Sqrt[7]}}

Pulling it all together, we get:

In[4]:= Graphics3D[{
           {Lighter@Lighter@Red,
            Sphere[centers, r]},

           {Lighter@Green,
            Sphere[{0, 0, 0}, R]},

           Opacity[0],
           EdgeForm[Thick],
           Parallelepiped[{-2, -2, -h/2}, 
            {{4, 0, 0}, {0, 4, 0}, {0, 0, h}}]} /. sol /. {r -> 1, R -> 2},

       Boxed -> False]

EDIT to add: I've updated this answer a couple times, each time to shift more work to Mathematica. In my first attempt, which I didn't post, I got the wrong answer because of a stupid typo when I plugged the radii of the spheres into the Pythagorean theorem; switching over to Norm fixed that. Then I told Solve to assume that $x$ and $z$ were positive (while using the "obvious" fact that $|y| = |x|$ to eliminate an equation before solving), and did some complicated thing with Tuples and Transpose to get all eight spheres. Eliminating those assumptions and recasting everything in terms of absolute values meant Solve took care of all of that automatically.

EDIT again to tweak the Solve expression to specify the domain for Reals and include y, which seems necessary for everything to work right in Mathematica v11.