  The path-following scheme in [Loisel and Maxwell, SIAM J. Matrix Anal. Appl.,
39-4 (2018), pp. 1726-1749] is adapted to efficiently calculate the dispersion
relation curve for linear surface waves on an arbitrary vertical shear current.
This is equivalent to solving the Rayleigh instability equation with linearised
free-surface boundary condition for each sought point on the curve. Taking
advantage of the analyticity of the dispersion relation, a path-following or
continuation approach is adopted. The problem is discretized using a
collocation scheme, parametrised along either a radial or angular path in the
wave vector plane, and differentiated to yield a system of ODEs. After an
initial eigenproblem solve using QZ decomposition, numerical integration
proceeds along the curve using linear solves as the Runge--Kutta $F(\cdot)$
function; thus, many QZ decompositions on a $2N$ companion matrix are exchanged
for one QZ decomposition and a small number of linear solves on a size $N$
matrix. A piecewise interpolant provides dense output. The integration
represents a nominal setup cost afterwhich very many points can be computed at
negligible cost whilst preserving high accuracy. Furthermore, a 2-dimensional
interpolant suitable for scattered data query points in the wave vector plane
is described. Finally, a comparison is made with existing numerical methods for
this problem, revealing that the path-following scheme is asymptotically two
orders of magnitude faster in number of query points.
