In this paper we investigate the stability of braneworld models constructed with non-minimally coupled multi-scalar fields. It is known that the tensor and vector perturbations are stable while the stability of the scalar perturbations are still unknown for such braneworld models. Models constructed with a single scalar are very different from those with multi-scalar fields. The linear scalar perturbed equations of the former can always be written as a supersymmetric Schrödinger equation, so they are stable at linear level. However, in general it is not true for the latter. With the aid of technics in cosmology and nodal theorem we present a systematic approach to deal with problems of the scalar perturbations in such braneworld models. It is find that a constraint equation arises and this constraint is not completely independent. The coupled Schrödinger-like equations are derived for the KK modes of the scalar perturbations. The stability and localization of the scalar zero modes are analyzed. In particular, we use these results to analyze the $f(R)$ braneworld model with an explicit solution and find that the scalar perturbations are stable and the scalar zero modes can not be localized on the brane, which ensure that there is no extra long-range force and the Newtonian potential on the brane can be recovered. ; Comment: 17 pages, 3 figures