We present and numerically implement a computational method to construct relativistic scattering amplitudes that obey analyticity, crossing, elastic and inelastic unitarity in three and four spacetime dimensions. The algorithm is based on the Mandelstam representation of the amplitude and iterations of unitarity. The input for the iterative procedure is given by the multi-particle double spectral density, the S-wave inelasticity, and the value of the amplitude at the crossing-symmetric point. The output, obtained at the fixed point of the iteration of unitarity, is a nonperturbative scattering amplitude. The amplitudes we obtain exhibit interesting features, such as non-zero particle production, intricate high-energy and near the two-particle threshold behavior. Scattering amplitudes obtained by initializing the iteration process with zero (or small) multi-particle input end up close to saturating the S-matrix bounds derived by other methods. There is a version of the iterative algorithm that is directly related to Feynman diagrams: it effectively re-sums infinitely many two-particle reducible planar Feynman graphs in the \phi^4 theory, which remarkably produces a unitary nonperturbative scattering amplitude function. Finally, we discuss how the algorithm can be further refined by including multi-particle unitarity.
