We discuss and implement a spectral method approach to computing stationary and axisymmetric black hole solutions and their properties in modified theories of gravity. The resulting code is written in the Julia language ...