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.