As in boiling super-heated liquids, the decay of a false vacuum is a first-order phase transition. A local ground state decays to an energetically more favorable minimum of lower energy due to the thermal and quantum fluctuations of the fields. In this work, we present an efficient semi-analytic method that computes the decay rate of such a state for any number of scalar fields and space-time dimensions. It is based on the collection of an arbitrary number of linear segments that describe a potential with several minima. The exact evolution of the field for each segment to provide the complete description of the bounce field configuration, which provides the leading contribution of the decay rate. By increasing the number of segments, one obtains the bounce action up to the desired precision. The resulting matching equations are solved semi-analytically and the generalization to more fields is computed iteratively via linear analytic perturbations. Based on this construction, we provide a robust and user-friendly Mathematica package that implements our method, named FindBounce. As it preserves the semianalytic structure of the method, its computational time grows linearly with the number of fields and segments. We present several applications and comparisons with other tools, where typical running time is roughly less than 1 (2) seconds for 10 (20) fields with 0.5% accuracy of the action. Finally, we describe a procedure that computes subleading contributions of the decay rate for any smooth potentials and extend it to include potentials with discontinuous first derivatives. As a consequence, we exhibit an exact decay rate at one loop for a real and complex scalar field in a bi-quartic potential with two treelevel minima. We compute the product of eigenvalues, remove the translational zero modes and renormalize the divergences with the zeta function formalism. We end up with a complete decay rate in a closed-form.