Solution of Boltzmann equation by a spherical-harmonic expansion approach is a computationally-efficient alternative to Monte Carlo. In this dissertation we extend this technique to compute the distribution function in multiple bands of silicon, using a multi-band band-structure which is accurate for high energies. A new variable transformation is applied on the spherical harmonic equations. This transformation (a) improves the numerical properties of the quations by enhancing the diagonal dominance of the resulting equations; (b) accounts for exponential dependence of the distribution function on energy as well as electric potential; and (c) opens the possibility of using superior Poisson solvers (d) while retaining the linearity of the original equations intact. The resulting Boltzmann equations are discretized using the current-conserving control-volume approach. The discretized equation are solved using line successive-over-relaxation (SOR) method. Numerical noise in the distribution was analyzed to be originating from the absence of coupling. Noise is removed by using acoustic phonons in inelastic approximation. A novel self-adjoint easy-to-discretize formulation for the inelastic acoustic phonons is developed. A test case of thermal equilibrium for multi-band is derived and used to validate the code. Hole-continuity and Poisson equation were solved along with the multi-band Boltzmann equations. The equations are solved in a Gummel-type decoupled loop. A \nnn\ device is simulated to test the simulator. The simulator is then applied to study a one-dimensional short-base bipolar junction transistor. While these simulations are self-consistent, a two-dimensional sub-micron MOSFET is simulated in a non-self-consistent manner.