An efficient algorithm is presented for the reliable computation of frequency response matrices of the form

. Efficiency is achieved through an initial reduction of

to upper Hessenberg form

so that (

) remains upper Hessenberg as ω is varied over a desired range. Linear systems with upper Hessenberg coefficient matrices can then be solved easily and accurately.