  We develop mathematical framework and computational tools for calculating
frequency responses of linear time-invariant PDEs in which an independent
spatial variable belongs to a compact interval. In conventional studies this
computation is done numerically using spatial discretization of differential
operators in the evolution equation. In this paper, we introduce an alternative
method that avoids the need for finite-dimensional approximation of the
underlying operators in the evolution model. This method recasts the frequency
response operator as a two point boundary value problem and uses
state-of-the-art automatic spectral collocation techniques for solving integral
representations of the resulting boundary value problems with accuracy
comparable to machine precision. Our approach has two advantages over currently
available schemes: first, it avoids numerical instabilities encountered in
systems with differential operators of high order and, second, it alleviates
difficulty in implementing boundary conditions. We provide examples from
Newtonian and viscoelastic fluid dynamics to illustrate utility of the proposed
method.
