  The detailed analysis of the global structure of the renormalization-group
(RG) flow diagram for a model with isotropic and cubic interactions is carried
out in the framework of the massive field theory directly in three dimensions
(3D) within an assumption of isotropic exchange. Perturbative expansions for RG
functions are calculated for arbitrary $N$ up to the four-loop order and
resummed by means of the generalized Pad$\acute{\rm e}$-Borel-Leroy technique.
Coordinates and stability matrix eigenvalues for the cubic fixed point are
found under the optimal value of the transformation parameter. Critical
dimensionality of the model is proved to be equal to $N_c=2.89 \pm 0.02$ that
agrees well with the estimate obtained on the basis of the five-loop
$\ve$-expansion [H. Kleinert and V. Schulte-Frohlinde, Phys. Lett. B342, 284
(1995)] resummed by the above method. As a consequence, the cubic fixed point
should be stable in 3D for $N\ge3$, and the critical exponents controlling
phase transitions in three-dimensional magnets should belong to the cubic
universality class. The critical behavior of the random Ising model being the
nontrivial particular case of the cubic model when N=0 is also investigated.
For all physical quantities of interest the most accurate numerical estimates
with their error bounds are obtained. The results achieved in the work are
discussed along with the predictions given by other theoretical approaches and
experimental data.
