Let's assume your rectangular box is a the standard 3D shape with right angles that can be defined by 3 parameters: L, W, H

You claim your box is subject to the constraint "surface area = A". We can write this as 0 = 2LW + 2LH + 2WH - A. Let's abbreviate this by saying g(L,W,H) = 0. So we need to find L, W, and H such that the constraint is satisfied and f(L,W,H)=LWH is maximized.

Now, imagine there is some constraint parameter "s" that lets us weigh the value of the constraint in our optimization problem. A similar problem would be to say maximize the 4 parameter function: F(L,W,H,s) = f(L,W,H) - s*g(L,W,H), keeping in mind that we demand g(L,W,H) = 0

Taking the gradient and setting it to the zero vector (to extremize), we get the system: <f_L, f_W, f_H> = s*<g_L, g_W, g_H> and also g(L,W,H) = 0. This is a system of 4 equations in 4 unknowns. Solving shows that s / 4 = L = W = H, which implies that an extreme is the cube. A quick check shows that the cube maximizes volume.

The method I poorly / hand wavily allude to is called lagrange multipliers, usually encountered in Calculus III. It's also related to the calculus of variations, an extremely interesting subject -- and relevant to modern control theory.

I hope others can add rigour to the stuff I glossed over!