  It was recently shown, that in a class of tensor-multi-scalar theories of
gravity with a nontrivial target space metric, there exist scalarized neutron
star solutions. An important property of these compact objects is that the
scalar charge is zero and therefore, the binary pulsar experiments can not
impose constraints based on the absence of scalar dipole radiation. Moreover,
the structure of the solutions is very complicated. For a fixed central energy
density up to three neutron star solutions can exist -- one general
relativistic and two scalarized, that is quite different from the scalarization
in other alternative theories of gravity. In the present paper we address the
stability of these solutions using two independent approaches -- solving the
linearized radial perturbation equations and performing nonlinear simulations
in spherical symmetry. The results show that the change of stability occurs at
the maximum mass models and all solutions before that point are stable. This
leads to the interesting consequence that there exists a stable part of the
scalarized branch close to the bifurcation point where the mass of the star
increases with the decrease of the central energy density.
