The Gaussian and Laguerre orthogonal ensembles are fundamental to random matrix theory, and the marginal eigenvalue distributions are basic observable quantities. Notwithstanding a long history, a formulation providing high precision numerical evaluations for $N$ large enough to probe asymptotic regimes, has not been provided. An exception is for the largest eigenvalue, where there is a formalism due to Chiani which uses a combination of the Pfaffian structure underlying the ensembles, and a recursive computation of the matrix elements. We augment this strategy by introducing a generating function for the conditioned gap probabilities. A finite Fourier series approach is then used to extract the sequence of marginal eigenvalue distributions as a linear combination of Pfaffians, with the latter then evaluated using an efficient numerical procedure available in the literature. Applications are given to illustrating various asymptotic formulas, local central limit theorems, and central limit theorems, as well as to probing finite size corrections. Further, our data indicates that the mean values of the marginal distributions interlace with the zeros of the Hermite polynomial (Gaussian ensemble) and a Laguerre polynomial (Laguerre ensemble).