  We study an inhomogeneous Neumann boundary value problem for functions of
least gradient on bounded domains in metric spaces that are equipped with a
doubling measure and support a Poincar\'e inequality. We show that solutions
exist under certain regularity assumptions on the domain, but are generally
nonunique. We also show that solutions can be taken to be differences of two
characteristic functions, and that they are regular up to the boundary when the
boundary is of positive mean curvature. By regular up to the boundary we mean
that if the boundary data is $1$ in a neighborhood of a point on the boundary
of the domain, then the solution is $-1$ in the intersection of the domain with
a possibly smaller neighborhood of that point. Finally, we consider the
stability of solutions with respect to boundary data.
