Polydispersity is believed to have important effects on the formation of liquid crystal phases in suspensions of rodlike particles. To understand such effects, we analyze the phase behavior of thin hard rods with length polydispersity. Our treatment is based on a simplified Onsager theory, obtained by truncating the series expansion of the angular dependence of the excluded volume. We describe the model and give the full phase equilibrium equations; these are then solved numerically using the moment free energy method which reduces the problem from one with an infinite number of conserved densities to one with a finite number of effective densities that are moments of the full density distribution. The method yields exactly the onset of nematic ordering. Beyond this, results are approximate but we show that they can be made essentially arbitrarily precise by adding adaptively chosen extra moments, while still avoiding the numerical complications of a direct solution of the full phase equilibrium conditions. We investigate in detail the phase behavior of systems with three different length distributions: a (unimodal) Schulz distribution, a bidisperse distribution, and a bimodal mixture of two Schulz distributions which interpolates between these two cases. A three-phase isotropic-nematic-nematic coexistence region is shown to exist for the bimodal and bidisperse length distributions if the ratio of long and short rod lengths is sufficiently large, but not for the unimodal one. We systematically explore the topology of the phase diagram as a function of the width of the length distribution and of the rod length ratio in the bidisperse and bimodal cases.