  This work proposes a set of equations that can be used to numerically compute
spacetimes containing a stationary black hole. The formalism is based on the
3+1 decomposition of General Relativity with maximal slicing and spatial
harmonic gauge. The presence of the black hole is enforced using the notion of
apparent horizon in equilibrium. This setting leads to the main result of this
paper: a set of boundary conditions describing the horizon and that must be
used when solving the 3+1 equations. Those conditions lead to a choice of
coordinates that is regular even on the horizon itself. The whole procedure is
validated with three different examples chosen to illustrate the great
versatility of the method. First, the single rotating black holes are recovered
up to very high values of the Kerr parameter. Second, non-rotating black holes
coupled to a real scalar field, in the presence of a negative cosmological
constant (the so-called MTZ black holes), are obtained. Last, black holes with
complex scalar hairs are computed. Eventually, prospects for future work, in
particular in contexts where stationarity is only approximate, are discussed.
