The existence of multi black hole solutions in asymptotically flat spacetimes is one of the expectation from the final state conjecture in general relativity. In this talk, I will present preliminary works in this direction via a semilinear toy model in dimension 3. In particular, I show 1) an algorithm to construct approximate solutions to the energy critical wave equation that converge to a sum of solitons at an arbitrary polynomial rate in (t-r); 2) a robust method to solve the remaining error terms for the nonlinear equation. The methods apply to energy supercritical problems. This work is part of my PhD thesis.