Numerical simulation of a three-dimensional semiconductor device is a fundamental problem in information science. The mathematical model is defined by a nonlinear system of initial-boundary problem including four partial differential equations: an elliptic equation for electrostatic potential, two convection-diffusion equations for electron concentration and hole concentration, a heat conduction equation for temperature. The electrostatic potential appears within the concentration equations and heat conduction equation, and the electric field strength controls the concentrations and the temperature. The electric field potential is solved by the conservative block-centered method, and the order of the accuracy is improved by the electric potential. The concentrations and temperature are computed by the upwind block-centered multistep method, where three different numerical methods are involved. The multistep method is adopted to approximate the time derivative. The block-centered method is used to discretize the diffusion. The upwind scheme is applied to approximate the convection to avoid numerical dispersion and nonphysical oscillation. The block-centered difference simulates diffusion, concentrations, temperature, and the adjoint vector functions simultaneously. It has the local conservation of mass, which is an important nature in numerical simulation of a semiconductor device. By using the variation, energy estimates, induction hypothesis, embedding theorem and the technique of a priori estimates of differential equations, convergence of the optimal order is obtained. Numerical examples are provided to show the effectiveness and viability. This method provides a powerful tool for solving the challenging benchmark problem.