In this paper, we study the braneworld scenarios in the presence of two real scalar fields coupled by gravity. The first-order formalism for the bent brane (for both de Sitter and anti-de Sitter geometry), leads us to discuss the shape invariance method in the bent brane systems. So, by using the fluctuations of metric and fields we obtain the Schr¨odinger equation. Then we factorize the corresponding Hamiltonian in terms of multiplication of the first-order differential operators. These first-order operators lead us to obtain the energy spectrum with the help of shape invariance method.
