In the present Letter, we study the braneworld scenarios in the presence of quintom dark energy 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 quintom fields we obtain the Schrödinger 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.