Two-dimensional incompressible Navier-Stokes equations are solved numerically to model the thermohydrodynamic performance of a dynamically loaded journal bearing which is modeled as eccentrically rotating cylinders. The region between those cylinders are occupied by Newtonian lubricants, whose physical properties such as viscosity and thermal conductivity are assumed to be functions of local temperature. A single domain pseudospectral method which combines Fourier expansions and Chebyshev polynomials for spatial discretization is introduced in conjunction with appropriate time marching scheme for the unsteady incompressible Navier-Stokes equations. The selection of these polynomial functions is favorable since both FFT algorithms for Fourier and Chebyshev expansions are easily available. In this numerical model, the journal is dynamically loaded by an external force and set free, so that its center moves in such a way to strike a balance between the applied load and the hydrodynamic forces. The pseudo-spectral scheme is then applied to a few classical problems, such as concentric rotating cylinders and journal bearings with lubricants of constant and varying (temperature dependent) viscosity to establish the validity of the numerical scheme in simulating these problems realistically as well as to gauge the convergence characteristics and relevant numerical issues. The numerical modeling has been found to be reasonably accurate and robust enough to serve as a tool for the study the flow in the region between the journal and the bearing.