recurrence relation, all terms of the sequence positive
Finally I resolved the issue of existence and uniqueness. We will regard each $a_n$ as rational function of $a_1$. That is, we will consider the sequence of rational functions
$$ a_0 = a_0(x) = 0, \qquad a_{1} = a_1(x) = x, \qquad a_{n+1} = \frac{n}{a_n} - a_n - a_{n-1}. \tag{*} $$
We would like to address the set of all $x\in\mathbb{R}$ for which $a_n(x) > 0$ for all $n \geq 1$. To this end, write
$$ I_n = \{ x \in \mathbb{R} : a_1(x) > 0, \cdots, a_n(x) > 0 \}, \qquad I_{\infty} = \bigcap_{n=1}^{\infty} I_n. $$
Then the solution set is exactly $I_{\infty}$, which we seek to identify. In this answer, we will prove that $I_{\infty}$ is a singleton, i.e.,
Proposition. There exists a unique $x \in \mathbb{R}$ for which $a_n(x) > 0$ for all $n \geq 1$.
We know that $I_1 = (0, \infty)$. Now we prove the following statement:
Claim. For any $n \geq 2$, the followings hold:
- $I_n = (\alpha, \beta)$ for some $0 \leq \alpha < \beta < \infty$. Moreover,
- If $n$ is even, then $a_{n-1}(\alpha) = 0$ and $a_n(\beta) = 0$,
- If $n$ is odd, then $a_{n-1}(\beta) = 0$ and $a_n(\alpha) = 0$.
- $(-1)^k a_{k+1}'(x) \geq (-1)^{k-1}a_k'(x) \geq 1$ for all $1 \leq k \leq n-1$ and $x \in I_n$.
Proof of Claim. We invoke induction.
(Base case) The claim is easy to verify directly when $n = 2$. Indeed, we know that $a_2(x) = \frac{1}{x}-x$, and so, $a_2(1) = 0$ and $I_2 = (0, 1)$. Also, $-a_2'(x) = \frac{1}{x^2} + 1 \geq 1$ and $a_1'(x) = 1$ proves part (2) of the claim for $n = 2$.
-
(Induction step) Assume that the claim holds for a given $n \geq 2$. Then by the induction hypothesis,
\begin{align*} (-1)^n a_{n+1}'(x) &= \left(\frac{n}{a_n^2} + 1\right) (-1)^{n-1}a_n'(x) - (-1)^{n-2}a_{n-1}'(x) \\ &\geq \frac{n}{a_n^2} (-1)^{n-1}a_n'(x) \\ & > 0. \end{align*}
Now we consider two case
$\text{Case 1.} \ $ If $n$ is even, then $a_n(\alpha^+) = +\infty$ and hence $a_{n+1}(\alpha^+) = -\infty$. Also we have $a_{n+1}(\beta^-) = +\infty$. Since $\alpha_{n+1}' > 0$ on $I_n$, there exists a unique zero $\gamma \in I_n$. Therefore $I_{n+1} = (\gamma, \beta)$.
$\text{Case 2.} \ $ If $n$ is odd, then by a similar argument, we can check that there exists a unique zero $\gamma$ of $a_{n+1}$ in $I_n$ and $I_{n+1} = (\alpha, \gamma)$.
So the part 1 of Claim for $n+1$ is proved. Finally, on $I_{n+1}$, we have $a_{n+1} > 0$, and so,
$$ \frac{n}{a_n^2} = 1 + \frac{a_{n+1} + a_{n-1}}{a_n} > 1. $$
Plugging this back shows that $(-1)^n a_{n+1}' \geq (-1)^{n-1}a_n'$ and therefore the induction step is proved. ////
Now write $I_n = (\alpha_n, \beta_n)$. Then by the intermediate step of the above proof, we can easily check that $(\alpha_n)_{n\geq 2}$ is non-decreasing, $(\beta_n)_{n\geq 2}$ is non-increasing and
$$ \alpha_{n} < \alpha_{n+2} < \beta_{n+2} < \beta_n. $$
So both $(\alpha_n)$ and $(\beta_n)$ converges. Moreover, if $\alpha = \lim \alpha_n$ and $\beta = \lim \beta_n$, then
$$ I_{\infty} = [\alpha, \beta]. $$
Finally, we prove that $\alpha = \beta$. To this end, we prove the following claim.
Claim. There exists $c \in (0, 1)$ such that $a_n(x) \leq c \sqrt{n}$ for all $n\geq 0$ and $x \in I_{\infty}$. In particular,
$$ |a_{n}'(x)| \geq \frac{1}{c^{2(n-1)}}. $$
Proof of Claim. Write $a_n = a_n(x)$ for simplicity. Then
$a_n^2 \leq n$ shows that $a_n \leq \sqrt{n}$ for all $n \geq 0$.
-
We have
$$\frac{n}{a_n} - a_n = a_{n-1} + a_{n+1} \leq \sqrt{n-1} + \sqrt{n+1} \leq 2\sqrt{n}. $$
Using this, we can easily conclude that $a_n \geq c_1\sqrt{n}$ for $c_1 = \sqrt{2}-1$.
-
Similarly, we find that
$$\frac{n}{a_n} - a_n = a_{n-1} + a_{n+1} \geq c_1(\sqrt{n-1} + \sqrt{n+1}) \geq c_1 \sqrt{2n}. $$
Here, we utilized the inequality $\sqrt{t-1}+\sqrt{t+1} \geq \sqrt{2t}$ which holds for all $t \geq 1$. Then by choosing $c$ as $c = \left(\sqrt{\smash[b]{2+c_1^2}} - c_1\right)/\sqrt{2} \approx 0.749$, we find that $a_n \leq c\sqrt{n}$ for all $n\geq 0$.
Finally, as in the proof of the previous claim, we note that
$$ |a_{n+1}'(x)| \geq \frac{n}{a_n^2} |a_n'(x)| \geq \frac{1}{c^2}|a_n'(x)|. $$
Therefore the desired claim follows. ////
Now we are ready to prove the uniqueness. Assume otherwise that $\alpha < \beta$. Then by the mean-value theorem,
$$|a_n(\alpha) - a_n(\beta)| = |a_n'(x)|(\beta - \alpha) \geq \frac{\beta-\alpha}{c^{2(n-1)}} $$
for some $x \in (\alpha, \beta)$. On the other hand, the left-hand side is bounded by
$$ |a_n(\alpha) - a_n(\beta)| \leq a_n(\alpha) + a_n(\beta) \leq 2c\sqrt{n}. $$
This is a contradiction as $n\to\infty$.
I post it as a curiosity. Considering the difference equation
$$ \frac{n}{x_n^2} = 1 + \frac{x_{n+1}+x_{n-1}}{x_n} $$
and considering
$$ \Delta x_n = x_n-x_{n-1} $$
we have
$$ \frac{n}{x_n^2} = 3+\frac{\Delta x_{n+1}-\Delta x_n}{x_n} $$
which is a difference approximation for the DE
$$ \frac{t}{x^2(t)} = 3 + \frac{\ddot x(t)}{x(t)} $$
now calling $x(1) = a = \frac{2 \Gamma \left(\frac{3}{4}\right)}{\Gamma \left(\frac{1}{4}\right)}$ and $x(2) = a-\frac 1a$ and integrating we have
In red the plot for $x_n = \sqrt{\frac n3 + \frac{1}{36n}+ \mathcal{O}(\frac{1}{n^3})}$ and in blue the DE solution.
and now considering instead the initial conditions $x(2) = a-\frac 1a, x(3) = \frac{1-3a^2}{a^3-3}$ we obtain
This result points in the direction of good agreement between the DE and the recurrence. Now following with the initial conditions
$$ x(3) = \frac{1-3a^2}{a^3-3}\\ x(4) = \frac{8 a^3-4 a}{3 a^4-4 a^2+1} $$
we obtain
The following answer is the special case $p=1$ of part 3 of this answer to a closely related question in which also the conjecture of Sangchul Lee's comment is proved.
Uniqueness: Suppose that $a_n,a_n'$ were two such sequences satisfying $a_0=a_0'=0$ and the recursion $$\label{eq3}\tag1 a_{n+1}=\frac n{a_n}-a_n-a_{n-1}$$ related to two different values $a_1=a,a_1'=a'.$ Then $d_n=a_n-a_n'$ satisfy $d_{n+1}=-\left(\frac n{a_na_n'}+1\right)d_n-d_{n-1}$ for $n\geq1$ and $d_0=0$. Now we have $a_na_{n+1}+a_n^2+a_na_{n-1}=n$ and hence $a_n\leq \sqrt n$. The same holds for $a_n'$: $a_n'\leq \sqrt n$. Hence $$|d_{n+1}|\geq 2|d_n|-|d_{n-1}|\mbox{ for }n\geq1.$$
We have $|d_2|\geq2|d_1|$ and now prove by induction that $$|d_{n}|\geq\frac n{n-1}|d_{n-1}|\mbox{ for }n\geq2.$$ This is true for $n=2$ and if it is true for some $n$ then $|d_{n-1}|\leq\frac{n-1}n |d_n|$ and hence $$|d_{n+1}|\geq 2|d_n|-|d_{n-1}|\geq\left(2-\frac{n-1}n\right)|d_n|=\frac{n+1}n|d_n|.$$
As a consequence, $|d_n|\geq n|d_1|=n|a-a'|$ for all $n$, which contradicts $|d_n|=|a_n-a_n'|\leq\sqrt n$. Therefore $a$ and $a'$ must be equal.
Existence: We use the function $w(n,z)$ defined for real $z$ and positive integer $n$ as the unique positive solution $w$ of $\frac nw-w=z$. It can be given explicitly as \begin{equation}\nonumber w(n,z)=-\frac z2+\sqrt{\frac{z^2}4+n}=\frac n{\frac z2+\sqrt{\frac{z^2}4+n}}. \end{equation} Observe that for any positive $n$, the mapping $z\to w(n,z)$ is strictly decreasing because the derivative of the mapping $w\to\frac nw-w$ is always negative.
We define the sequence $U_1(n)=\sqrt{n}$ for all $n$. Then we define recursively sequences $L_k,U_k$ by $L_k(0)=U_k(0)=0$ for all $k$ and \begin{equation}\nonumber L_k(n)=w(n,U_k(n+1)+U_k(n-1))\mbox{ and }U_{k+1}(n)=w(n,L_k(n+1)+L_k(n-1)) \end{equation} for all $n,k\geq1.$
By induction and using that $z\mapsto w(n,z)$ is strictly decreasing it follows that $$L_k(n)\leq L_{k+1}(n)\leq U_{k+1}(n)\leq U_k(n)$$ for all $n,k$.
Thus the pointwise limits
\begin{equation}\nonumber
U(n)=\lim_{k\to\infty}U_k(n)\mbox{ and } L(n)=\lim_{k\to\infty}L_k(n)
\end{equation}
exist because for fixed $n$, the sequences $U_k(n)$ and $L_k(n)$, $k=1,2,3,...$ are monotonous and
bounded. The properties of the sequences $U_k,L_k$ imply besides $U(0)=L(0)=0$
a) $L(n)\leq U(n)$ for every $n$,
b) $U(n)=w(n,L(n+1)+L(n-1))$ and $L(n)=w(n,U(n+1)+U(n-1))$ for all positive $n$.
As a consequence of b), the two sequences $A_n,B_n$, $n=0,1,...$, defined by
$A_n=U(n),\,B_n=L(n)$ if $n$ is even and
$A_n=L(n),\,B_n=U(n)$ if $n$ is odd
both satisfy the recursion (\ref{eq3}). They must be equal because we have proved uniqueness of positive solution sequences of (\ref{eq3}).
Observe that the above construction of $U,L$ can be used to approximate the positive sequence numerically.