on the powerful method of the superpotential [47–49], which can be useful even for two or more scalar fields (see, for example, [50] in the context of soft-wall models). However, even if one solution is constructed in this way, there is no guarantee that this is the only solution with the same action. We will now describe how to look for all possible solutions of a given action using a combination of numerical and graphical techniques. Whenever there is more than one static solution to the above boundary value problem with the same action , we say that multiple solutions exist. In general, the bulk scalar fields can have Dirichlet boundary conditions, Neumann boundary conditions or more general mixed boundary conditions. Here we focus on the case where we have one kink field φ (obeying Dirichlet boundary conditions), with the remaining N − 1 fields χ$_{a}$ having Neumann or mixed boundary conditions. When the profiles of these extra fields are monotonic, they will tend to stabilize the extra dimension, whereas if their profiles have vanishing derivatives inside the interval, they will tend to destabilize the extra dimension [42]. Despite this subtlety, we will generically refer to the non-kink fields as “stabilization” fields. To find solutions we proceed as follows: we specify the Lagrangian in the bulk and on one of the branes, and we numerically solve an initial value problem to determine the profiles of the fields along the extra dimension. Dirichlet boundary conditions are imposed on the kink field φ at the initial brane by demanding that it vanish there. For this to hold, we assume the kink field has a sufficiently heavy brane mass so that it decouples from the stabilization fields on the branes. As a result, the kink field disappears from the junction conditions (35)-(36), which then yield only N conditions on the initial brane. This leaves N + 1 initial conditions that need to be specified, which we take to be the boundary values for the derivatives φ $^{′}$, χ ′ 1 . . . χ ′ N − $_{1}$, and σ $^{′}$. After solving the initial value problem for a given choice of initial conditions, we impose Dirichlet boundary conditions on φ at the final boundary by locating the second brane at a point where the profile of φ vanishes. In general, the profile will vanish at several points along the extra dimension, and one may study kinks with the desired number of nodes by choosing the location of the second brane accordingly. Here, as in the flat case, we are primarily interested in nodeless kink solutions, and we therefore place the second brane at the first zero of the profile function. 14 A. Multiple Solutions