  We present a systematically improvable method for numerically solving
relativistic three-body integral equations for the partial-wave projected
amplitudes. The method consists of a discretization procedure in momentum
space, which approximates the continuum problem with a matrix equation. It is
solved for different matrix sizes, and in the end, an extrapolation is employed
to restore the continuum limit. Our technique is tested by solving a three-body
problem of scalar particles with an $S$ wave two-body bound state. We discuss
two methods of incorporating the pole contribution in the integral equations,
both of them leading to agreement with previous results obtained using
finite-volume spectra of the same theory. We provide an analytic and numerical
estimate of the systematic errors. Although we focus on kinematics below the
three-particle threshold, we provide numerical evidence that the methods
presented allow for determination of amplitude above this threshold as well.
