  In this work we study the stability regions of linear multistep or
multiderivative multistep methods for initial-value problems by using
techniques that are straightforward to implement in modern computer algebra
systems. In many applications, one is interested in (i) checking whether a
given subset of the complex plane (e.g. a sector, disk, or parabola) is
included in the stability region of the numerical method, (ii) finding the
largest subset of a certain shape contained in the stability region of a given
method, or (iii) finding the numerical method in a parametric family of
multistep methods whose stability region contains the largest subset of a given
shape. First we describe a simple procedure to exactly calculate the stability
angle $\alpha$ in the definition of $A(\alpha)$-stability. As an illustration,
we consider two finite families of implicit multistep methods: we exactly
compute the stability angles for the $k$-step BDF methods ($3\le k\le 6$) and
for the $k$-step second-derivative multistep methods of Enright ($3\le k\le
7$). Next we determine the exact value of the stability radius in the BDF
family for each $3\le k\le 6$, that is, the radius of the largest disk in the
left half of the complex plane, symmetric with respect to the real axis,
touching the imaginary axis and lying in the stability region of the
corresponding method. Finally, we demonstrate how some Schur--Cohn-type
theorems of recursive nature and not relying on the RLC method can be used to
exactly solve some optimization problems within infinite parametric families of
multistep methods. As an example, we choose a two-parameter family of
implicit-explicit (IMEX) methods: we identify the unique method having the
largest stability angle in the family, then we find the unique method in the
same family whose stability region contains the largest parabola.
