  In this paper, we show that the eigenvalues and eigenvectors of the spectral
discretisation matrices resulted from the Legendre dual-Petrov-Galerkin (LDPG)
method for the $m$th-order initial value problem (IVP): $u^{(m)}(t)=\sigma
u(t),\, t\in (-1,1)$ with constant $\sigma\not=0$ and usual initial conditions
at $t=-1,$ are associated with the generalised Bessel polynomials (GBPs). The
essential idea of the analysis is to properly construct the basis functions for
the solution and its dual spaces so that the matrix of the $m$th derivative is
an identity matrix, and the mass matrix is then identical or approximately
equals to the Jacobi matrix of the three-term recurrence of GBPs with specific
integer parameters. This allows us to characterise the eigenvalue distributions
and identify the eigenvectors. As a by-product, we are able to answer some open
questions related to the very limited known results on the collocation method
at Legendre points (studied in 1980s) for the first-order IVP, by reformulating
it into a Petrov-Galerkin formulation. Moreover, we present two stable
algorithms for computing zeros of the GBPs, and develop a general space-time
spectral method for evolutionary PDEs using either the matrix diagonalisation,
which is restricted to a small number of unknowns in time due to the
ill-conditioning but is fully parallel, or the QZ decomposition which is
numerically stable for a large number of unknowns in time but involves
sequential computations. We provide ample numerical results to demonstrate the
high accuracy and robustness of the space-time spectral methods for some
interesting examples of linear and nonlinear wave problems.
