  A code for the numerical evaluation of hyperelliptic theta-functions is
presented. Characteristic quantities of the underlying Riemann surface such as
its periods are determined with the help of spectral methods. The code is
optimized for solutions of the Ernst equation where the branch points of the
Riemann surface are parameterized by the physical coordinates. An exploration
of the whole parameter space of the solution is thus only possible with an
efficient code. The use of spectral approximations allows for an efficient
calculation of all quantities in the solution with high precision. The case of
almost degenerate Riemann surfaces is addressed. Tests of the numerics using
identities for periods on the Riemann surface and integral identities for the
Ernst potential and its derivatives are performed. It is shown that an accuracy
of the order of machine precision can be achieved. These accurate solutions are
used to provide boundary conditions for a code which solves the axisymmetric
stationary Einstein equations. The resulting solution agrees with the
theta-functional solution to very high precision.
