The steady-state Reynolds equation for gas lubricating films leads to a quasilinear boundary-value problem in two dimensions. This equation contains a first-order derivative term whose coefficient is nonlinear and also large for most of the practical cases. Thus, this becomes a quasilinear singular perturbation problem. Existing numerical schemes are faced with the failure of convergence as well as numerical instability which often results in overshooting the answers. Asymptotic approximations also give poor results if the film-thickness ratios are not small. In this paper an accurate and reliable numerical scheme is presented. Convergence is proved independent of the bearing numbers and film-thickness ratios. A weighted upwind-difference form is used to discretize the differential equation. Theory of M-matrices and associated inequalities are employed to prove the ensuing monotonicity. The analysis presented in this paper can be extended to a more general class of singular-perturbation quasilinear boundary-value problems. Numerical results and graphs for pressure distributions and bearing loads are provided for the parabolic slider. Comparisons are made with other existing results concerning numerical as well as asymptotic analysis.