We present a numerical method to solve the neutrino Boltzmann equation coupled to stellar core collapse and lepton conservation, with no approximations made to the neutrino scattering kernels. Spherical symmetry is assumed. Our finite differencing of the Boltzmann equation is similar to its finite difference representation in the discrete ordinates method, but our auxiliary equations relating the zone-center and zone-edge distribution functions are different. We also differ from the discrete ordinates method in the way we solve the finite differenced Boltzmann equation. With our particular choice of auxiliary relations, we avoid oscillations in the neutrino distribution function that are characteristic of other choices, and we also obtain the correct diffusion limit without the introduction of ad hoc terms in the radial transport. Our finite differenced Boltzmann equation is fully implicit with regard to time centering. It is in conservative form, and upwind differencing is used for the radial and angular transport. Gaussian quadratures are used for the neutrino direction cosines. The scattering terms are explicitly written in terms of the neutrino distribution function, which is being solved for. They are not assumed known. The finite differenced Boltzmann equation, the finite differenced electron fraction equation (the electron fraction equation is derived from the equation for lepton conservation), and the finite differenced equation that describes the change in the specific internal energy of the fluid from absorption and emission of electron neutrinos, which is operator split from the specific internal energy equation, are linearized and solved using a multidimensional Newton-Raphson iteration. The linearized equations take on a block tridiagonal form and are solved using the Feautrier elimination scheme. Because we are concerned with radiation transport in hydrodynamic media, we also present methods for finite differencing the angular aberration and frquency shift terms in the Boltzmann equation. In particular we present a finite difference representation for the angular aberration term suitable for use with Gaussian quadrature. We also present a generalization of Bruenn's method for the frequency shift, suitable for use with a Boltzmann equation rather than a multigroup flux-limited diffusion approximation. This method is conservative for both neutrino number and neutrino energy. The matter in the stellar core to which the neutrinos are coupled is evolved using Newtonian gravity and 0(v/c) Lagrangian hydrodynamics. The finite differencing for the hydrodynamics equations is also presented. All of our numerical methods are designed to accommodate realistic equations of state.