  We propose a novel numerical method for solving multi-dimensional, special
relativistic Boltzmann equations for neutrinos coupled to hydrodynamics
equations. It is meant to be applied to simulations of core-collapse
supernovae. We handle special relativity in a non-conventional way, taking
account of all orders of v/c. Consistent treatment of advection and collision
terms in the Boltzmann equations is the source of difficulties, which we
overcome by employing two different energy grids: Lagrangian remapped and
laboratory fixed grids. We conduct a series of basic tests and perform a
one-dimensional simulation of core-collapse, bounce and shock-stall for a
15M_{sun} progenitor model with a minimum but essential set of microphysics. We
demonstrate in the latter simulation that our new code is capable of handling
all phases in core-collapse supernova. For comparison, a non-relativistic
simulation is also conducted with the same code, and we show that they produce
qualitatively wrong results in neutrino transfer. Finally, we discuss a
possible incorporation of general relativistic effects in our method.
