So let us first define what is a bound state, in particular with respect to scattering theory. Scattering theory studies asymptotic states coming out of interactions in the limit of infinite times. That means that if we want the bound state to appear in our out states, it must be stable, or not spontaneously decaying after a finite amount of time. If it does decay after a finite amount of time, it cannot be a part of an out state. Any other state is a resonance in the scattering matrix and can be, in principle, treated by perturbation theory.
The reason why bound states are not discussed in your typical QFT course is that 1) they do not occur in the models discussed in introductory courses ($\phi^4$ theory and QED), 2) there are very few "nice" pedagogical computations of bound states.
However, be it an "approximate bound state" (a resonance) or a proper bound state of 2 particles (the case of more particle is easy to generalize), one can find its properties by considering the Bethe-Salpeter equation (the wiki is not very well written but I link it nonetheless). The Bethe-Salpeter equation was derived by Salpeter and Bethe in 1951 (with some precursors by Gell-Mann and Low) from the requirement that the normalization of the total momentum of the two-particle system is below the sum of the masses the individual particles. One then obtains the following formal integral equation of the two-point wave function $\psi(x_1,x_2)$ whose stationary solutions represent bound states and resonances:
$$\psi(x_1,x_2) = -i \int G^{\rm free}_1(x_1,x'_1)G^{\rm free}_2(x_1,x'_2)V(x'_1,x'_2,y_1,y_2) \psi(y_1,y_2) d^4x'_1 d^4x'_2 d^4y_1 d^4y_2$$
where I have suppressed space-time indices $x^\mu$, and where $G^{\rm free}$ are free-field propagators and $V$ an interaction kernel consisting of a sum of all 4-point irreducible Feynman diagrams. The energies of the stationary solutions represent the energies of the bound states/resonances. By solving this equation iteratively and including only tree-level diagrams, you will get the "non-relativistic"/weak-field limit equivalent to a Schrödinger type problem in the first iteration. From higher iterations you can get QFT corrections to these types of problems.
If you are interested in the formation of these states, the wave-function actually also represents the transition amplitude of two free particles into a bound state, $\psi(x_1,x_2) = \langle \Psi_{\rm B}|T[\phi(x_1),\phi(x_2)]|\Omega_{\rm vac} \rangle$, which you can use to add an effective interaction vertex with a possible free outgoing bound state in your S-matrix computations. (Note, however, that the bound state necessarily has a smaller total momentum than the two particles, so a third particle such as a photon must carry energy away in the processes.)
