  The complete analysis of a model with three quartic coupling constants
associated with an O(2N)--symmetric, a cubic, and a tetragonal interactions is
carried out within the three-loop approximation of the renormalization-group
(RG) approach in $D=4-2\epsilon$ dimensions. Perturbation expansions for RG
functions are calculated using dimensional regularization and the minimal
subtraction (MS) scheme. It is shown that for $N\ge 2$ the model does possess a
stable fixed point in three dimensional space of coupling constants, in
accordance with predictions made earlier on the base of the lower-order
approximations. Numerical estimate for critical (marginal) value of the order
parameter dimensionality $N_c$ is given using Pad\'e-Borel summation of the
corresponding $\epsilon$--expansion series obtained. It is observed that
two-fold degeneracy of the eigenvalue exponents in the one-loop approximation
for the unique stable fixed point leads to the substantial decrease of the
accuracy expected within three loops and may cause powers of $\sqrt{\epsilon}$
to appear in the expansions. The critical exponents $\gamma$ and $\eta$ are
calculated for all fixed points up to $\epsilon^3$ and $\epsilon^4$,
respectively, and processed by the Borel summation method modified with a
conformal mapping. For the unique stable fixed point the magnetic
susceptibility exponent $\gamma$ for N=2 is found to differ in third order in
$\epsilon$ from that of an O(4)--symmetric point. Qualitative comparison of the
results given by $\epsilon$--expansion, three-dimensional RG analysis,
non-perturbative RG arguments, and experimental data is performed.
