NEW (Oct 27, 2017)
I can now show that the recurrence given in the Question is correct.
In fact, I can deduce a simpler and shorter recurrence relation.
For this, we consider the Kummer Surface $K$ associated to $J$.
This is a quartic surface in ${\mathbb P}^3$, and there is a map
$\kappa \colon J \to K$
that identifies a point $P$ and its negative $-P$ (so it is 2-to-1 except
at the 16 points of order 2 on $J$). In terms of the coordinates one
usually uses, the image of your point on $K$ is $(0:1:0:0)$. Also, the
constant term in the first component of the Mumford representation of a
point with image $(x_1:x_2:x_3:x_4)$ on $K$ is $x_3/x_1$ (assuming
$x_1 \neq 0$). Now there is a $4 \times 4$ matrix $B(\vec{x},\vec{y})$
of bi-quadratic forms in two sets of variables $\vec{x} = (x_1,\ldots,x_4)$
and $\vec{y} = (y_1,\ldots,y_4)$ (i.e., separately homogeneous of degree 2
in each set of variables) such that, if $\vec{x},\vec{y},\vec{w},\vec{z}$
are coordinate vectors of the images of points $P, Q, P-Q, P+Q$, resp.,
then, up to a common scaling factor,
$$B_{ii}(\vec{x},\vec{y}) = w_i z_i \quad\text{and}\quad
B_{ij}(\vec{x},\vec{y}) = w_i z_j + w_j z_i \quad\text{for $i \neq j$.}$$
(See the book by Cassels and Flynn for background.)
We write $\kappa(nP) = (1 : x_{n,2} : x_{n,3} : x_{n,4})$
(assuming the first coordinate is non-zero), so that $x_{n,3}$ is the term
of interest. It is easy to see that the matrix $B(\vec{x},\vec{y})$
determines the unordered pair $\{\vec{w}, \vec{z}\}$ (up to scaling
of the coordinate vectors).
Now we use the relations
$B(\vec{x}_n, \vec{x}_1) \sim \vec{x}_{n-1} * \vec{x}_{n+1}$
(where $\vec{w} * \vec{z}$ denotes the symmetric matrix with entries
$w_i z_i$ on the diagonal and $w_i z_j + w_j z_i$ off the diagonal)
for two successive values of $n$, together with the equations
$f(\vec{x}_{n-1}) = \ldots = f(\vec{x}_{n+2}) = 0$, where $f = 0$
is the equation defining $K$. Together they define an ideal in the
polynomial ring in 12 variables over $\mathbb Q$,
and we can compute the elimination ideal $I$ with
respect to $\{x_{n-1,3}, x_{n,3}, x_{n+1,3}, x_{n+2,3}\}$.
Magma says that $I$ is generated by 10 polynomials, one of which is
linear in $x_{n+2,3}$ (and also in $x_{n-1,3}$). Writing $a_n$ for
$x_{n,3}$ and solving, we find the three-term recurrence
$$a_{n+2} = a_n \frac{a_{n-1}+a_{n+1}}{a_{n-1} - a_{n+1}}
- 2 \frac{a_n^2 a_{n+1}^2 + a_n a_{n+1} + 5000(a_n + a_{n+1}) - 1}
{(a_{n-1}-a_{n+1}) a_n^2 a_{n+1}^2}
$$
We can also check that the four-term recurrence given in the question
(in terms of the algebraic relation it implies between five successive
values of $a_{n}$) is implied by the relations that hold between
five successive coordinate vectors.
In the spirit of a previous remark (that I removed when re-writing the
answer to give a complete solution), note that the variety $V_k$ given by
$$B(\vec{x}_n,\vec{x}_1) \sim \vec{x}_{n-1} * \vec{x}_{n+1}, \ldots,
B(\vec{x}_{n+k},\vec{x}_1) \sim \vec{x}_{n+k-1} * \vec{x}_{n+k+1}$$
and
$$f(\vec{x}_{n-1}) = \ldots = f(\vec{x}_{n+k+1}) = 0$$
is two-dimensional for any $k \ge 0$. So we can expect a relation between
any three quantities depending on the variables. For example, there is a
relation relating $a_{n-1}$, $a_n$ and $a_{n+1}$, but it is of degree 4
in $a_{n-1}$ and in $a_{n+1}$ and of degree 6 in $a_n$. We need to consider
more quantities to find a relation that is linear in one of them.
C.jacobian(). $\endgroup$