Numerical solvers of the incompressible Navier-Stokes equations have reproduced turbulence phenomena such as the law of the wall, the dependence of turbulence intensities on the Reynolds number, and experimentally observed properties of turbulence energy production. In this article, we begin a sequence of investigations whose eventual aim is to derive and implement numerical solvers that can reach higher Reynolds numbers than is currently possible. Every time step of a Navier-Stokes solver in effect solves a linear boundary value problem. The use of Green's functions leads to numerical solvers which are highly accurate in resolving the boundary layer, which is a source of delicate but exceedingly important physical effects at high Reynolds numbers. The use of Green's functions brings with it a need for careful quadrature rules and a reconsideration of time steppers. We derive and implement Green's function based solvers for the channel flow and plane Couette flow geometries. The solvers are validated by reproducing turbulence phenomena in good agreement with earlier simulations and experiment.