Non-linear rational recurrence related to the Jacobian of hyperelliptic curve EDIT The initial revision had typo, hope it is fixed.
Let $C$ be the hyperelliptic curve $y^2=x^5+x+5^2$ and $J$ its Jacobian.
Let $P=J(0,5)$ and define $a(n)$ to be the constant coefficient of
the first coordinate of $nP$.
$a(n)$ starts $1, 0, 0, 1/5000, 12500, 200000003/125000$.
Experimentally $a(n)$ satisfy the nonlinear recurrence up to $n=100$:
a(n + 2) = -((50000002*a(n - 1)*a(n - 2) - 10000*a(n - 2) - 5)*a(n + 1)*a(n) - 5*(2000*a(n - 2) + 1)*a(n - 1)*a(n) - 5*a(n - 1)*a(n - 2) + 25000*a(n - 2) + 100000014)/(10000*(5000*a(n - 2) - 1)*a(n - 1)*a(n) + (50000002*a(n - 1)*a(n - 2) + 2*(25000001*a(n - 1) + 25000000*a(n - 2) - 5000)*a(n) - 10000*a(n - 2) - 5)*a(n + 1) - 10000*a(n - 1)*a(n - 2) + 25000)

Is the recurrence true?
sagemath code:
def jacob1(KK,lim=100):
    """
    K=QQ
    """
    pr=PolynomialRing(KK,'x')
    x=pr.gen()
    w0=KK(5)
    dg=5
    pol= x**5+x+5**2  
    C = HyperellipticCurve(pol)
    pt0=C([KK(0),KK(5)**1,KK(1)])
    J=C.jacobian()
    P=J(KK)(pt0)
    f=[]
    for i in srange(lim):
        Q= i*P
        u=Q[0]
        T=u.constant_coefficient()
        f.append(T)
        print u
        print '  ',T
    return f

 A: 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.
