The idea is to use the node Distribute Points on Faces to generate points in a cube volume, the reason for using this node is that it has an option for minimum distance when in Poisson-disk mode (I believe it does that by generating the points and then deleting the invalid ones), which we can use to keep spheres from intersecting.
Since the node Distribute Points on Faces only distributes on faces, we can simulate a volume with many layers of planes inside the volume. With a bigger layer count, more possible positions are available to points, and thus points tend to be more closer to each other, but always obeying the minimum distance.
to eliminate the need to check and delete spheres intersecting with the bounding cube, the layers can be created with edges far from the cube bounds by the desired radius, that way spheres will at most touch the bounding cube.
First a vertical line is created centered to the bounding cube, it's length is $w - 2r$, where $w$ is the bounding cube's size, and $r$ is the desired radius for the spheres.
Then, planes with dimensions equal to $w - r$ are instanced equally spaced along the line. Spheres generated on them will be more closer to each other as the number of planes increase.
On the generated planes, points are distributed with a minimum distance of $2r$. (points of all planes are considered when generating them, so the minimum distance also applies between points of different layers).
After that just instance spheres with radius $r$ on the distributed points;
I also check if a sphere of the desired radius can fit in the bounding cube before outputting the geometry, since a single sphere is still created by this method when the radius value meet that condition.