Integral representations for the solutions of the Laplace and modified Helmholtz equations can be obtained using Green's theorem. However, these representations involve both the solution and its normal derivative on the boundary, and for a well-posed boundary-value problem (BVP) one of these functions is unknown. Determining the Neumann data from the Dirichlet data is known as constructing the Dirichlet-to-Neumann map. A new transform method was introduced in Fokas (1997, Proc. R. Soc. Lond. A, 53, 1411–1443) for solving BVPs for linear and integrable nonlinear partial differential equations (PDEs). For linear PDEs this method can be considered as the analogue of the Green's function approach in the Fourier plane. In this method the Dirichlet-to-Neumann map is characterized by a certain equation, the so-called global relation, which is formulated in the complex k-plane, where k denotes the complex extension of the spectral (Fourier) variable. Here we solve the global relation numerically for the Laplace and modified Helmholtz equations in a convex polygon. This is achieved by evaluating the global relation at a properly chosen set of points in the spectral (Fourier) plane, which is why this method has been called a ‘spectral collocation method’. Numerical experiments suggest that the method inherits the order of convergence of the basis used to expand the unknown functions, namely, exponential for a polynomial basis such as Chebyshev, and algebraic for a Fourier basis. However, the condition number of the associated linear system is much higher for a polynomial basis than for a Fourier one.