This paper considers a Monte-Carlo Nystrom method for solving integral equations of the second kind, whereby the values (z(yi))(1 <= i <= N) of the solution z at a set of N random and independent points (y(i))1(<i<N) are approximated by the solution (zN,i)(1 <= i <= N) of a discrete N dimensional linear system obtained by replacing the integral with the empirical average over the samples (yi)(1 <= i <= N). Under the unique assumption that the integral equation admits a unique solution z(y), we prove the invertibility of the linear system for sufficiently large N with probability one, and the convergence of the solution (zN,i)(1 <= i <= N) towards the point values (z(yi))(1 <= i <= N) in a mean-square sense at a rate O(N-1/2 ). For a particular choices of kernels, the discrete linear system arises as the Foldy-Lax approximation for the scattered field generated by a system of N sources emitting waves at the points (yi)(1 <= i <= N). In this context, our result can be equivalently considered as a proof of the well-posedness of the Foldy-Lax approximation for systems of N point scatterers, and of its convergence as N -> +infinity in a mean-square sense to the solution of a Lippmann-Schwinger equation characterizing the effective medium. The convergence of Monte-Carlo solutions at the rate O(N-1/2) is numerically illustrated on one-dimensional examples and for solving a two-dimensional Lippmann-Schwinger equation.