Based upon a variational principle and the associated theory derived in three preceding papers, an expression for the magneto-elastic buckling value for a system of an arbitrary number of parallel superconducting beams is given. The total current is supposed to be equal both in magnitude and direction for all beams, and the cross-sections are circular. The expression for the buckling value is formulated more explicitly in terms of the so-called buckling amplitudes, the latter following from an algebraic eigenvalue problem. The pertinent matrix is formulated in terms of complex functions, which are replaced by real potentials. The matrix elements are calculated by a numerical method, solving a set of integral equations with regular kernels. Apart from the buckling value(s) the buckling modes are also obtained. Finally, our results are compared with the results of a mathematically less complicated theory, i.e. the method of Blot and Savart.