The numerical properties of a deterministic Boltzmann equation solver based on a spherical harmonics expansion of the distribution function are analyzed and improved. A fully coupled discretization scheme of the Boltzmann and Poisson equations is proposed, where stable equations are obtained based on the H-transformation. It is explicitly shown that the resultant Jacobian matrix for the zeroth order component has property M for a first order expansion, which improves the stability even of higher order expansions. The detailed dependence of the free-streaming operator and the scattering operator on the electrostatic potential is exactly considered in the Newton-Raphson scheme. Therefore, convergence enhancement is achieved compared with previous Gummel-type approaches. This scheme is readily applicable to small-signal and noise analysis. As numerical examples, simulation results are shown for a silicon n + nn + structure including a magnetic field, an SOI NMOSFET and a SiGe HBT.