This system seems simple enough to throw to a CAS or some optimizer for small $n$. The inequalities are sequential: once you have understood behaviour for $n$, you can add a few more constraints for $n+1$ and try to understand the expanded system. To illustrate, I use $x_n$ in place of $f(n)$.

The first round is to maximize $x_0$ subject to $x_0^2 + x_0 \leq 1$. This is solved by the original poster and results in $x_0 = 1/\phi$.

The next round ($n=2$) adds the constraint equivalent to $x_{n-1} \leq 1/(1+x_0)$, and a little calculus and checking of boundary conditions gives $x_1 + x_0$ is maximized at $x_1=x_0=1/\phi$.

With each increase of $n$, the additional constraints are of the form $x_n \leq 1-x_j x_{n-j}$. However, we can ignore those constraints and focus on maximizing $x_0 + x_n$, which again we get with $x_0=x_n=1/\phi$. So I am guessing $1/\phi$ is the best you can do for the average value. (This does not quite work. However, if you get a value early on that approaches $1/\phi$, then this impacts later values. More work needs to be done.)

Here is some more work. Let us suppose we have $n$ least that $x_n \gt 1/\phi$. Then we have all the $x_i \it 1/\phi$ for all $0 \leq i \lt n$, and we know by above that the sum of the first $n+1$ variables is at most $(n+1)/\phi$. In particular, if $x_n=1/\phi + \epsilon$, then $x_0 \lt 1/\phi \lt \epsilon$.

If we now consider $x_j = 1/\phi - \delta$ and $x_{n-j} = 1/\phi - \gamma$, the relevant inequality gives $\phi\epsilon \lt \gamma + \delta$, so for $n \gt 2$ we have the sum is less than $(n+1)/\phi$ by at least $(\phi -1)\epsilon$, and this is magnified for all $\lfloor n/2 \rfloor$ pairs. It remains to study the case when $x_1$ and perhaps when $x_2$ are greater than $1/\phi$.

**Update 2018.12.28**

This argument is pretty simple, and suggests that something similar might work for $f$ over lower bounded subsets of the rationals or reals.

One sees that constant $f$ with value $f \leq 1/\phi$ satisfies the system of inequalities, so the best we can hope for is that $d_n \leq c=1/\phi$. I will recoordinatize and set $e_n = f(n) - c$, and use $c$ to simplify the TeX. We have $cc + c = 1$ and can rewrite the functional relation as $e_j(c + e_{n-j})+ ce_{n-j} + e_n \leq 0$.

We will be interested in showing the sum of the $e_j$ is at most 0. Towards that we rewrite the inequality as $e_{n-j} + e_n \leq e_j(-c - e_{n-j}) + (1 - c)e_{n-j}$. If all the $e_k$ are nonpositive, we are done. Otherwise, pick $j$ least with $e_j \gt 0$, and let $n\geq j$ be least $n$ with $e_n \gt 0$. Then the first form of the equality will imply $e_{n+j} \lt 0$, and the second form shows $-c \leq e_{n-j} \lt e_{n-j} + e_n \lt 0$.

Now let $j$ be least with $e_j \gt 0$. For every $n$, either $e_n \leq 0$ or else $e_{n-j} + e_n\lt 0$. This implies that the finite sum $e_0 + e_1 + \cdots + e_n \leq 0$ for all choices satisfying all the inequalities.

Gerhard "That's How I See It" Paseman, 2018.12.22.