[Edited to add sharper bound (number of sign changes) and connection
with "Descartes' Rule of Signs"]
Yes, the zeros at $x=1,2,\ldots,i-1$ are the only real zeros of $f_i$.
We prove that in general an "exponential polynomial" with $d+1$ nonzero
terms, i.e. $A(x) = \sum_{j=1}^{d+1} a_j \exp(\lambda_j x)$
with distinct $\lambda_j \in {\bf R}$ and each $a_j \in {\bf R}^*$,
can have at most $d$ real roots, counted with multiplicity.
We use induction on $d$, the base case $d=0$ being trivial.
Suppose we've proved the case $d-1$ for some $d>0$.
Now $A_0(x) := e^{-\lambda_1 x} A(x)$ is an exponential polynomial,
with the same number of nonzero terms and the same roots as $A$,
whose $j=1$ term is constant. Hence $A'_0(x) := \frac{d}{dx}(A_0(x))$ is an
exponential polynomial with only $d$ nonzero terms. By the inductive
hypothesis, it thus has at most $d-1$ real roots. But by Rolle's
theorem (easily generalized to allow multiple roots), there's at least
one root of $A'_0(x)$ between each consecutive pair of roots of $A_0$, and
thus of $A$. Therefore $A$ can have no more than $d$ roots with
multiplicity. This completes the induction step and the proof.
In the case at hand, $d=i-1$, each $\lambda_j = \log j$, and each $a_j
= (-1)^{i-j} {i \choose j}$. Having already located $i-1=d$ distinct
real roots of $A$, we deduce that there are no others (and that each of
the known roots is simple). QED
[added later] It is no accident that in our case, where $A(x)$
actually has $d$ distinct roots, the coefficients $a_j$ alternate
in sign when the $\lambda_j$ are listed in increasing order.
Indeed if we assume without loss of generality that
$\lambda_1 < \lambda_2 < \cdots < \lambda_{d+1}$
then the number of real roots (counted with multiplicity) is at most
the number, call it $s$, of sign changes in the coefficient sequence
$(a_1,a_2,\ldots,a_{d+1})$; that is, at most the number of
$j \in \lbrace 1,2,\ldots,d \rbrace$ such that $a_j a_{j+1} < 0$.
To see this, we argue by induction as above, noting also that
$A_0(x) = \sum_{j=1}^{d+1} a_j \exp((\lambda_j - \lambda_1) x)$
has $\lambda_j - \lambda_1 > 0$, so each coefficient
$(\lambda_j - \lambda_1) a_j$ of $A'_0(x)$ has the same sign as $a_j$
for $j>1$. Thus if $a_1 a_2 < 0$ then $A'_0(x)$ has $s-1$ sign changes,
and we're done as before. If $a_1 a_2 > 0$ then $A'_0(x)$
has $s$ sign changes. But then $A_0$ is monotone on $x \leq x_0$
where $x_0$ is the smallest root of $A'_0(x)$.
Since the two leading terms of $A_0$ as $x \rightarrow -\infty$
have the same sign, it follows that $A_0$ does not vanish on
$x \leq x_0$, and thus by Rolle that $A_0$ has no more real zeros
than its derivative. Therefore in this case too $A_0$ has at most $s$
real zeros. This completes the inductive step and the proof. QED
This argument may feel familiar, most likely because the same technique
proves Descartes' "rule of signs", which bounds the number of positive
roots of an ordinary real polynomial by its number of sign changes.
Indeed this "rule of signs" is equivalent to the special case
$\lambda_j \in \lbrace 0, 1, 2, \ldots \rbrace$ of our bound,
when $A(x)$ is an ordinary polynomial in $e^x$.
Here is a solution to your warm up problem. It uses a few known elementary identities, and some short inductions for identities I didn't recognize. I also changed your notation slightly from $l$ to $k$ in the internal summations.
Setting $x=2$, the first term vanishes, and we combine the second and third terms as $(B)$ and take the last term as $(A)$.
\begin{align}
&\sum_{j=0}^i \binom{n-1-j}{i-j} 2^{j-1} \sum_{k=0}^j (-1)^k \binom{n+1}{k} 2^{i-k-1}(i-k) \tag{A} \\
&\sum_{j=0}^{i-2} \binom{n-1-j}{i-2-j} 2^{j-1} \sum_{k=0}^j (-1)^k \binom{n+1}{k} \left( 2^{i-k} + 2^{i-k-1}(i-2-k)\right) \tag{B}
\end{align}
Fortunately for us, both are equal to $i\cdot2^{i-2}$ when $i$ is even, and do not depend on $n$.
First, we focus on $(A)$. Switch the order of summation to get
\begin{equation}
\sum_{k=0}^i (-1)^k \binom{n+1}{k} 2^{i-2}(i-k) \sum_{j=k}^i \binom{n-1-j}{i-j} 2^{j-k} \tag{A}
\end{equation}
Then reindex the internal summation, and repeatedly apply the hockey stick identity to show
\begin{equation}
\sum_{j=0}^{i-k}\binom{n-1-k-j}{i-k-j}2^j = \sum_{j=0}^{i-k}\binom{n-k}{j}
\end{equation}
I found it easier to write out by changing variables $a=i-k-j$, $b=i-k$ and $c=n-i-1$, then simplifying $\sum_{a=0}^b \binom{c+a}{c}2^{b-a}$.
$(A)$ naturally splits at the point $i-k$, and since we want to show that the whole thing is $i\cdot2^{i-2}$, we can break it into two pieces and show
\begin{align}
\sum_{k=1}^i (-1)^k k\binom{n+1}{k}\sum_{j=0}^{i-k} \binom{n-k}{j} &=0 \tag{A1}\\
\sum_{k=0}^i (-1)^k \binom{n+1}{k} \sum_{j=0}^{i-k} \binom{n-k}{j} &=1 \tag{A2}
\end{align}
We can simplify $(A1)$ slightly by incorporating $k$ into the binomial, and ignoring the $n+1$ term that comes out.
Switching the order of summations again, we write
\begin{equation}
Q(n,i) = \sum_{j=0}^{i-1} \sum_{k=1}^{i-j} (-1)^k \binom{n}{k-1} \binom{n-k}{j}
\end{equation}
We will induct on $n$ and $i$ to show that $Q(n,i)=0$ for $i$ even, and $-1$ for $i$ odd. The base cases are easy to check, especially if taking $i=0$ and $i=1$.
But first, we need an auxillary identity, which we will also use in $(A2)$.
\begin{equation}
P(n,i) = \sum_{j=0}^{i} (-1)^{i-j}\binom{n}{j} \binom{n-1-j}{i-j}
\end{equation}
We claim that $P(n,i)=1$ for all $n$ and $i$. This is certainly true for $i=0$. We split $\binom{n}{j}$ into two to get an induction:
\begin{align}
P(n,i) &= \sum_{j=0}^{i} (-1)^{i-j}\binom{n-1}{j} \binom{n-1-j}{i-j} + \sum_{j=1}^{i} (-1)^{i-j}\binom{n-1}{j-1} \binom{n-1-j}{i-j} \\
&= \sum_{j=0}^{i} (-1)^{i-j}\binom{n-1}{i} \binom{i}{j} + \sum_{j=1}^{i} (-1)^{i-j}\binom{n-1}{j-1} \binom{n-1-j}{i-j} \\
&= \binom{n-1}{i}\sum_{j=0}^{i} (-1)^{i-j}\binom{i}{j} + \sum_{j=0}^{i-1} (-1)^{i-1-j}\binom{n-1}{j} \binom{(n-1)-1-j}{(i-1)-j} \\
&= 0 + P(n-1,i-1)
\end{align}
Now we give a similar proof for $Q$. In the fourth to fifth lines, we used the alternating sum of binomial coefficients up to some number.
You might be able to give a direct proof if your generatingfunctionology is strong, since these look like convolutions of simple functions, but these proofs seemed easy enough that I didn't bother trying.
\begin{align}
Q(n,i) &= \sum_{j=0}^{i-1} \sum_{k=1}^{i-j} (-1)^k \binom{n}{k-1} \binom{n-k}{j} \\
&= \sum_{j=0}^{i-1} \sum_{k=1}^{i-j} (-1)^k \binom{n-1}{k-1} \binom{n-k}{j} + \sum_{j=0}^{i-2} \sum_{k=2}^{i-j} (-1)^k \binom{n-1}{k-2} \binom{n-k}{j} \\
&= \sum_{j=0}^{i-1} \sum_{k=1}^{i-j} (-1)^k \binom{n-1}{j} \binom{n-1-j}{k-1} + \sum_{j=0}^{i-2} \sum_{k=2}^{i-j} (-1)^k \binom{n-1}{k-2} \binom{n-k}{j} \\
&= \sum_{j=0}^{i-1} \binom{n-1}{j} \sum_{k=0}^{i-1-j} (-1)^{k+1} \binom{n-1-j}{k} + \sum_{j=0}^{i-2} \sum_{k=1}^{i-1-j} (-1)^{k+1} \binom{n-1}{k-1} \binom{n-1-k}{j} \\
&= \sum_{j=0}^{i-1} \binom{n-1}{j} (-1)^{i-1-j+1} \binom{n-2-j}{i-1-j} - Q(n-1, i-1) \\
&= -P(n-1,i-1) - Q(n-1, i-1)
\end{align}
Now, since our equation $(A1)$ was just $(n+1)Q(n,i)$, and $i$ is even, it's $0$. A similar treatment yields $(A2)$, again using the alternating binomial sum identity near the end:
\begin{align}
&\,\sum_{k=0}^i (-1)^k \binom{n+1}{k} \sum_{j=0}^{i-k} \binom{n-k}{j} \\
&= \sum_{j=0}^i \sum_{k=0}^{i-j} (-1)^k \binom{n+1}{k} \binom{n-k}{j} \\
&= \sum_{j=0}^i \sum_{k=0}^{i-j} (-1)^k \binom{n}{k} \binom{n-k}{j} + \sum_{j=0}^i \sum_{k=0}^{i-j} (-1)^k \binom{n}{k-1} \binom{n-k}{j} \\
&= \sum_{j=0}^i \sum_{k=0}^{i-j} (-1)^k \binom{n}{j} \binom{n-j}{k} + Q(n,i) \\
&= \sum_{j=0}^i \binom{n}{j} \sum_{k=0}^{i-j} (-1)^k \binom{n-j}{k} \\
&= \sum_{j=0}^i \binom{n}{j} (-1)^{i-j} \binom{n-1-j}{i-j} \\
&= P(n,i)
\end{align}
So we have shown that $(A) = i\cdot2^{i-2}$ for $i$ even. Let's use this for $(B)$, since $i-2$ is also even, we have:
\begin{equation}
\sum_{j=0}^{i-2} \binom{n-1-j}{i-2-j} 2^{j-1} \sum_{k=0}^j (-1)^k \binom{n+1}{k} 2^{i-k-3}(i-2-k) = (i-2)2^{i-4}
\end{equation}
After multiplying both sides by $2^2$, the only place the left side differs from $(B)$ is the extra $-2$ in $(i-2-k)$, and the right side is $2^{i-1}$ smaller than we would like.
So we show that
\begin{equation}
\sum_{j=0}^{i-2} \binom{n-1-j}{i-2-j} 2^{j-1} \sum_{k=0}^j (-1)^k \binom{n+1}{k} 2^{i-k} = 2^{i-1}
\end{equation}
Dividing out the factor of $2^{i-1}$ and switching the order of summation, we get
\begin{equation}
\sum_{k=0}^{i-2} (-1)^k \binom{n+1}{k} \sum_{j=k}^{i-2} \binom{n-1-j}{i-2-j}2^{j-k}
\end{equation}
Of course we recognize our hockey stick identity from earlier, so this simplifies to the case of $(A2)$
\begin{equation}
\sum_{k=0}^{i-2} (-1)^k \binom{n+1}{k} \sum_{j=0}^{i-2-k} \binom{n-k}{j} =1
\end{equation}
Best Answer
Everything is already contained in OEIS comments for A001864 and A000435 (a remarkable comment is that A000435 is the sequence that started it all: the first sequence in the database!)
We take $n$ labelled vertices, consider all trees on them, and sum up the distances between all pairs of vertices (each distance counted twice).
One way to do it is the following: this sum is the number of 5-tuples $(T,a,b,c,d)$ such that $T$ is a tree, $a,b,c,d$ are vertices, $ab$ is an edge of $T$ and this edge belongs to the path between $c$ and $d$ (in the order $cabd$ on the path). If we remove $ab$, we get two connected components $A\ni a$, $B\ni b$. If $|A|=i$, $|B|=n-i$, we may fix $A$, $B$ by $\binom{n}i$ ways, after that fix restrictions of $T$ onto $A$, $B$ by $i^{i-2}(n-i)^{n-i-2}$ ways and fix $a,b,c,d$ by $i^2(n-i)^2$ ways. Totally we get RHS of your formula.
Why we get LHS is explained in Claude Lenormand's comment for A000435 (there we count the sum of distances from the fixed vertex 0 to other vertices in all trees, of course it is $n$ times less than the sum of all distances.)