The existence of multi black hole solutions in General Relativity is one of the expectations from the final state conjecture, the analogue of soliton resolution. 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.