The integral equation for the flow velocity u(x;k) in the steady Couette flow derived from the linearized Bhatnagar–Gross–Krook–Welander kinetic equation is studied in detail both theoretically and numerically in a wide range of the Knudsen number k between 0.003 and 100.0. First, it is shown that the integral equation is a Fredholm equation of the second kind in which the norm of the compact integral operator is less than 1 on Lp for any 1≤p≤∞ and thus there exists a unique solution to the integral equation via the Neumann series. Second, it is shown that the solution is logarithmically singular at the endpoints. More precisely, if x=0 is an endpoint, then the solution can be expanded as a double power series of the form ∑n=0∞∑m=0∞cn,mxn(xlnx)m about x=0 on a small interval x∈(0,a) for some a>0. And third, a high-order adaptive numerical algorithm is designed to compute the solution numerically to high precision. The solutions for the flow velocity u(x;k), the stress Pxy(k), and the half-channel mass flow rate Q(k) are obtained in a wide range of the Knudsen number 0.003≤k≤100.0; and these solutions are accurate for at least twelve significant digits or better, thus they can be used as benchmark solutions.