  We present a computational approach to the classical Schottky problem based
on Fay's trisecant identity for genus $g\geq 4$. For a given Riemann matrix
$\mathbb{B}\in\mathbb{H}^{g}$, the Fay identity establishes linear dependence
of secants in the Kummer variety if and only if the Riemann matrix corresponds
to a Jacobian variety as shown by Krichever. The theta functions in terms of
which these secants are expressed depend on the Abel maps of four arbitrary
points on a Riemann surface. However, there is no concept of an Abel map for
general $\mathbb{B} \in \mathbb{H}^{g}$. To establish linear dependence of the
secants, four components of the vectors entering the theta functions can be
chosen freely. The remaining components are determined by a Newton iteration to
minimize the residual of the Fay identity. Krichever's theorem assures that if
this residual vanishes within the finite numerical precision for a generic
choice of input data, then the Riemann matrix is with this numerical precision
the period matrix of a Riemann surface. The algorithm is compared in genus 4
for some examples to the Schottky-Igusa modular form, known to give the Jacobi
locus in this case. It is shown that the same residuals are achieved by the
Schottky-Igusa form and the approach based on the Fay identity in this case. In
genera 5, 6 and 7, we discuss known examples of Riemann matrices and
perturbations thereof for which the Fay identity is not satisfied.
