This paper describes a new method to calculate the stationary scattering matrix and its derivatives for Euclidean waveguides. This is an adaptation and extension to a procedure developed by Levitin and Strohmaier which was used to compute the stationary scattering matrix on surfaces with hyperbolic cusps (Levitin and Strohmaier, 2019), but limited to those surfaces. At the time of writing, these procedures are the first and only means to explicitly compute such objects. In this context the challenge we faced was that on Euclidean waveguides, the scattering matrix naturally inhabits a Riemann surface with a countably infinite number of sheets making it more complicated to define and compute. We overcame this by breaking up the waveguide into compact and non-compact components, systematically describing the resolvent for the Neumann Laplace operator on both of them, giving a thorough treatment of the Riemann surface, and then using a “gluing” construction (Melrose, 1995) to define the resolvent on the whole surface. From the resolvent, we were able to obtain the scattering matrix. The algorithm we have developed to do this not only computes the scattering matrix itself on such domains, but also arbitrarily high derivatives of it directly. We have applied this, together with the finite element method, to calculate resonances for a selection of domains and will present the results of some numerical calculations in the final section. Whilst this is certainly not the first, nor only method to compute resonances on these domains, i.e. Levitin and Marletta have done so previously (Levitin and Marletta, 2008) and Aslanyan, Parnovski and Valiliev before them (Aslanyan et al., 2000) and other techniques, such as perfectly matched layers may be adapted for this purpose (Jiang and Xiang, 2020). The method described here has several advantages in terms of speed and accuracy and moreover, provides more information about the scattering phenomena.