We present a method of calculating the interacting S-matrix to an arbitrary perturbative order for a large class of boson interaction Lagrangians. The method takes advantage of a previously unexplored link between the n-point Green's function and a certain system of linear Diophantine equations. By finding all nonnegative solutions of the system, the task of perturbatively expanding an interacting S-matrix becomes elementary for any number of interacting fields, to an arbitrary perturbative order (irrespective of whether it makes physical sense) and for a large class of scalar boson theories. The method does not rely on the position-based Feynman diagrams and promises to be extended to many perturbative models typically studied in quantum field theory. Aside from interaction field calculations we showcase our approach by expanding a pair of Unruh–DeWitt detectors coupled to Minkowski vacuum to an arbitrary perturbative order in the coupling constant. We also link our result to Hafnian as introduced by Caianiello and present a method to list all (2n−1)!! perfect matchings of a complete graph on 2n vertices.