In this work, we compute the Wigner distribution function on one-dimensional devices from wave functions generated by solving the Schrödinger equation. Our goal is to investigate certain issues that we encountered in implementing Wigner transport equation solvers, such as the large discrepancies observed between the boundary conditions and the solution in the neighborhood of the boundaries. By evaluating the Wigner function without solving the Wigner transport equation, we intend to ensure that the actual boundary conditions are consistent with those commonly applied in literature. We study both single-and double-barrier unbiased structures. We use simple potential profiles, so that we can compute the wave functions analytically for better accuracy. We vary a number of structure geometry, material, meshing, and numerical parameters, among which are the contact length, the barrier height, the number of incident wave functions, and the numerical precision used for the computations, and we observe how the Wigner function at the device boundaries is affected. For the double-barrier structures, we look at the density matrix function and we study a model for the device transmission spectrum which helps explain the lobelike artifacts that we observe on the Wigner function.