  The spectral density for random matrix $\beta$ ensembles can be written in
terms of the average of the absolute value of the characteristic polynomial
raised to the power of $\beta$, which for even $\beta$ is a polynomial of
degree $\beta(N-1)$. In the cases of the classical Gaussian, Laguerre, and
Jacobi weights, we show that this polynomial, and moreover the spectral density
itself, can be characterised as the solution of a linear differential equation
of degree $\beta+1$. This equation, and its companion for the resolvent, are
given explicitly for $\beta=2$ and $4$ for all three classical cases, and also
for $\beta=6$ in the Gaussian case. Known dualities for the spectral moments
relating $\beta$ to $4/\beta$ then imply corresponding differential equations
in the case $\beta=1$, and for the Gaussian ensemble, the case $\beta=2/3$. We
apply the differential equations to give a systematic derivation of recurrences
satisfied by the spectral moments and by the coefficients of their $1/N$
expansions, along with first-order differential equations for the coefficients
of the $1/N$ expansions of the corresponding resolvents. We also present the
form of the differential equations when scaled at the hard or soft edges.
