Mathematically, the quantity you're asking about is called the operator norm. Unfortunately, there's no simple formula for it. If it's a fully general affine transformation - for instance, if it could have an arbitrary combination of rotations and nonuniform scales, in any order - then I'm afraid there's nothing for it but to use singular value decomposition. If you apply SVD to your matrix then the largest singular value will be the maximum radius of the resulting ellipsoid. The other singular values will also be its other two radii, and the SVD procedure can also extract the orientation of the axes for you.
Implementing SVD is not for the faint of heart, as it involves finding eigenvalues. If all you want is the singular values themselves, they are the square roots of the eigenvalues of M^T * M. So if you have a 3x3 eigenvalue solver handy, or you don't mind writing one, you can use that. If you want to extract the orientations of the axes as well, then it gets more involved as you have to find eigenvectors too. On that Wikipedia article there is a list of links to libraries for doing SVD, one of which you may be able to use in your project.
If the form of your matrix is restricted in such a way that nonuniform scale happens at most once and is the first transform applied, i.e. is rightmost when you're using column vectors, then you can simplify this to just look at the lengths of the transformed axis vectors. In that case alone - i.e. a single nonuniform scale followed by any sequence of rotations, reflections, and uniform scales - looking at just the axis vectors will give you the right answer.