  Detailed observations of phenomena involving black holes, be it via
gravitational waves or more traditional electromagnetic means, can probe the
strong field regime of the gravitational interaction. The prediction of
features in such observations requires detailed knowledge of the black hole
spacetime, both within and outside of General Relativity. We present here a new
numerical code that can be used to obtain stationary solutions that describe
black hole spacetimes in a wide class of modified theories of gravity. The code
makes use of a relaxed Newton-Raphson method to solve the discretized field
equations with a Newton's polynomial finite difference scheme. We test and
validate this code by considering static and spherically symmetric black holes
both in General Relativity, as well as in scalar-Gauss-Bonnet gravity with a
linear (linear scalar-Gauss-Bonnet) and an exponential
(Einstein-dilaton-Gauss-Bonnet) coupling. As a by-product of the latter, we
find that analytic solutions obtained in the small coupling approximation are
in excellent agreement with our fully non-linear solutions when using a linear
coupling. As expected, differences arise when using an exponential coupling. We
then use these numerical solutions to construct a fitted analytical model,
which we then use to calculate physical observables such as the innermost
stable circular orbit and photon sphere and compare them to the numerical
results. This code lays the foundation for more detailed calculations of black
hole observables that can be compared with data in the future.
