  We extend a recently developed numerical code to obtain stationary,
axisymmetric solutions that describe rotating black hole spacetimes in a wide
class of modified theories of gravity. The code utilizes a relaxed
Newton-Raphson method to solve the full nonlinear modified Einstein's Equations
on a two-dimensional grid with a Newton polynomial finite difference scheme. We
validate this code by considering static and axisymmetric black holes in
General Relativity. We obtain rotating black hole solutions in
scalar-Gauss-Bonnet gravity with a linear (linear scalar-Gauss-Bonnet) and an
exponential (Einstein-dilaton-Gauss-Bonnet) coupling and compare them to
analytical and numerical perturbative solutions. From these numerical
solutions, we construct a fitted analytical model and study observable
properties calculated from the numerical results.
