  This manuscript deals with the analysis of numerical methods for the full
discretization (in time and space) of the linear heat equation with Neumann
boundary conditions, and it provides the reader with error estimates that are
uniform in time. First, we consider the homogeneous equation with homogeneous
Neumann boundary conditions over a finite interval. Using finite differences in
space and the Euler method in time, we prove that our method is of order 1 in
space, uniformly in time, under a classical CFL condition, and despite its lack
of consistency at the boundaries. Second, we consider the nonhomogeneous
equation with nonhomogeneous Neumann boundary conditions over a finite
interval. Using a tailored similar scheme, we prove that our method is also of
order 1 in space, uniformly in time, under a classical CFL condition. We
indicate how this numerical method allows for a new way to compute steady
states of such equations when they exist. We conclude by several numerical
experiments to illustrate the sharpness and relevance of our theoretical
results, as well as to examine situations that do not meet the hypotheses of
our theoretical results, and to illustrate how our results extend to higher
dimensions.
