  Via a generalization of the pseudospectral method for numerical solution of
differential equations, a family of nonlinear algebraic identities satisfied by
the zeros of a wide class of orthogonal polynomials is derived. The
generalization is based on a modification of pseudospectral matrix
representations of linear differential operators proposed in the paper, which
allows these representations to depend on two, rather than one, sets of
interpolation nodes. The identities hold for every polynomial family
$\{p_\nu(x)\}_{\nu=0}^\infty$ orthogonal with respect to a measure supported on
the real line that satisfies some standard assumptions, as long as the
polynomials in the family satisfy differential equations $\mathcal{A} p_\nu(x)
=q_\nu(x) p_\nu(x)$, where $\mathcal{A}$ is a linear differential operator and
each $q_\nu(x)$ is a polynomial of degree at most $n_0 \in \mathbb{N}$; $n_0$
does not depend on $\nu$. The proposed identities generalize known identities
for classical and Krall orthogonal polynomials, to the case of the nonclassical
orthogonal polynomials that belong to the class described above. The
generalized pseudospectral representations of the differential operator
$\mathcal{A}$ for the case of the Sonin-Markov orthogonal polynomials, also
known as generalized Hermite polynomials, are presented. The general result is
illustrated by new algebraic relations satisfied by the zeros of the
Sonin-Markov polynomials.
