Numerical Investigation of the Stable Nocturnal Boundary Layer The governing equations for the wind field and temperature field within the flat nocturnal atmospheric boundary layer are a highly nonlinear system of parabolic PDEs. This system is discretized into a crude two-layer numerical model via the finite difference approximation and the Monin-Obukhov similarity theory, and analyzed as a set of ODEs. The steady state problem is also transformed into an equivalent system of first order ODEs and then discretized into a very accurate "multi-layer" model using the orthogonal collocation method. Some numerical techniques for nonlinear problems such as numerical continuation and bifurcation analysis are used to study the steady state solutions as some physical parameters vary. The resulting bifurcation diagrams from the two-layer and multi-layer models have qualitatively similar behavior. This implies that the two-layer model, though mathematically crude, does capture some essential features of the original system. Time dependent solutions of the two layer model are computed via the fourth-order Runge-Kutta technique, for various combinations of parameters, and they match and support related bifurcation diagrams. Physically realistic wind and temperature profiles over the boundary layer are computed from the "multi-layer" model. Our results imply that operational application of this type of model of frost or pollution dispersion may not be made with confidence for certain parameter regimes, and they have important implications for the predictability of the nocturnal boundary layer for frost prediction or pollution dispersion. Space discretization for simple parabolic PDEs from an AUTO demo via pseudo-spectral method with Chebyshev basis functions is very accurate, and seems promising for future application to our problem.