Although a solution has been posted using analytical means, I note that it is possible to see the nonnegativity of the diagonal elements of $B - B^2$ purely algebraically, via combinatorial interpretation of the entries of the matrix $B$.
Fix $n$ and regard $A$ as a function of its off-diagonal entries $a_{ij}$, $i \neq j$, with its diagonal entries determined hy the requirement that the row sums of $A$ are all zero. A monomial in the variables $a_{ij}$, $i \neq j$, can be thought of as a directed graph on the vertices $\{1,2,\dots,n\}$ without "self-loops": the power of $a_{ij}$ in the monomial records how many directed edges go from $i$ to $j$ in the graph (the distinctness of the indices $i$ and $j$ in the variables $a_{ij}$ means that no vertex is ever connected to itself by a single edge). When $G$ is such a directed graph I will write $a_G = \prod_{(i,j) \in G} a_{ij}$ for the corresponding monomial (this product being interpreted "with multiplicity" whenever $G$ has more than one edge from $i$ to $j$ - although such $G$ will not arise in the cases of interest below). A forest on the vertices $\{1,2,\dots,n\}$ is an undirected graph on these vertices whose connected components contain no loops (i.e., are trees), and a rooted forest is a forest whose connected components each have a distinguished vertex. A rooted forest can be regarded as a directed graph by orienting each edge toward the root of its component. I will write $\mathcal{F}_k$ for the set of all rooted forests on $\{1,\dots,n\}$ having exactly $k$ components (equivalently, having exactly $n-k$ edges), and for $1 \leq j \leq n$ I will write $\mathcal{F}_k(j)$ for the set of rooted forests on $\{1,\dots,n\}$ in which $j$ is a root. Note that there is exactly one element $F$ of $\mathcal{F}_n$ (the graph having vertices $\{1,2,\dots,n\}$ and no edges), and in this case the corresponding monomial $a_F$ is $1$ and the sets $\mathcal{F}_n(j)$ are all equal to $\mathcal{F}_n$.
Letting $B := x(xI - A)^{-1}$ and denoting the row $i$, column $j$ entry of $B$ by $B_{ij}$, and letting $d(xI - A) := \det(xI - A)/x$ (note that it is clear a priori that $x$ divides $\det(xI - A)$), the combinatorial interpretations of the entries of $B$ promised up above are that
$$
d(xI-A) = \sum_{k=1}^n \left(\sum_{F \in \mathcal{F}_k} a_F\right) x^{k-1},
$$
and that
$$
d(xI-A) B_{ij} = \sum_{k=1}^n \left(\sum_{F \in \mathcal{F}_k(j), i \sim_F j} a_F\right) x^{k-1}, \qquad 1 \leq i, j \leq n,
$$
where $i \sim_F j$ means that $i$ and $j$ belong to the same connected component of $F$ (a condition that is vacuous when $i=j$). In this last formula note that when $i \neq j$ the inside sum corresponding to the index value $k=n$ is empty (because the one element of $\mathcal{F}_n$ has no edges) and does not contribute to the sum.
As the OP observed, it is straightforward to show that $B$ is (right) stochastic, i.e., that it has all its row sums equal to $1$. The formulas above are a refinement of that observation. That the rows of $B$ sum to $1$, or equivalently that the rows of $d(xI-A) B$ sum to $d(xI - A)$, reflects the fact that for any row index $i$, the rooted forests on the vertices $\{1,\dots,n\}$ can be partitioned according to the root of the component that contains $i$, with $d(xI - A) B_{ij}$ accounting for those forests in which $i$ lies in a component in which $j$ is a root, and with $d(xI-A)$ accounting for all of these forests.
To see how these formulas imply the nonnegativity of the diagonal of $B - B^2$, note that for any stochastic matrix $B$ (i.e. not just $B = x(xI - A)^{-1}$ with $A$ as above) we have
$$
B_{ii} = \left(\sum_{j=1}^n B_{ij}\right) B_{ii} = \sum_{j=1}^n B_{ij} (B_{ji} + B_{ii} - B_{ji}) = [B^2]_{ii} + \sum_{j \neq i} B_{ij} (B_{ii} - B_{ji}).
$$
For $B$ of the above specific type (but not for all stochastic matrices), our formulas will show that $B_{ii} - B_{ji} \geq 0$ for all $j \neq i$, which then (by the formula displayed immediately above) implies $B_{ii} \geq [B^2]_{ii}$ for all $i$.
Indeed, while $d(xI - A) B_{ii}$ is a sum over all rooted forests in which $i$ is a root, for $i \neq j$ the term $d(xI-A) B_{ji}$ accounts for only those forests in which $i$ is a root and $j$ is in the same component as $i$. For $i \neq j$ the difference $p_{ji} := d(xI - A) B_{ii} - d(xI - A) B_{ji}$ is thus the sum
$$
p_{ji} = \sum_{k=1}^n \left(\sum_{F \in \mathcal{F}_k(i), j \not \sim_F i} a_F\right) x^{k-1};
$$
which is plainly nonnegative whenever $x$ and the $a_{ij}$ are because it is a sum of monomials in these variables with nonnegative integer coefficients (indeed, any positive coefficient is $1$). Since $d(xI - A)$ is also such a sum (and is positive for positive $x$ and positive values of the $a_{ij}$), the nonnegativity of $B_{ii} - B_{ij}$ for $i \neq j$, and hence the nonnegativity of $B_{ii} - [B^2]_{ii}$ for all $1 \leq i \leq n$ follows.
The formula for $d(xI-A) = \det(xI-A)/x$ and the formulas for the matrix entries $B_{ij}$ (which are themselves determinant formulas by e.g. Cramer's rule) are slightly nontrivial to verify. Note, for example, that while it is straightforward to see that $d(xI - A)$ is a sum of monomials in $x$ and $a_{ij}$ of degree $n-1$, it is perhaps less clear which monomials of this form occur in the expansion and which do not. When $n=4$, for example, the formula immediately shows that the monomial $a_{12} a_{24} a_{41}$ does not appear with nonzero coefficient in $\det(xI - A)/x$ (the corresponding directed graph has a loop), although the monomial $a_{12} a_{24} a_{41} x$ does indeed appear with nonzero coefficient in the product of the diagonal entries of $xI - A$ (as do all degree $4$ products of $x$ with three "$a$" variables taken from distinct rows). So the formula above for $d(xI-A)$ implies that a substantial amount of cancellation would take place were one to simply evaluate the determinant via Laplace expansion (for example).
Some general remarks on the formula for $d(xI - A)$:
- It is a polynomial of degree $n-1$ in $x$ with leading term $x^{n-1}$.
- The coefficient of $x^{n-2}$ in $d(xI-A)$ is the sum of all variables $a_{ij}$, $i \neq j$.
- The coefficient of $x^{n-3}$ in $d(xI-A)$ is the sum of all products of pairs of variables $a_{ij} a_{kl}$, $i \neq j$, $k \neq l$, where we additionally require that $a_{ij}$ and $a_{kl}$ not occupy positions that are reflections of one another across the diagonal, or in the same row - or more suggestively for the general case, where the directed edges from $i$ to $j$ and from $k$ to $l$ do not together form a loop (preventing the monomial from representing part of a tree) or point in two "directions" from the same point (preventing the monomial from representing part of a rooted tree).
- In general, the coefficient of $x^{n-(d+1)}$ in $d(xI-A)$ is a sum of degree $d$ monomials $a_{i_1 j_1} a_{i_2 j_2} \cdots a_{i_d j_d}$ in the variables $a_{ij}$, where a given monomial appears if and only if the directed graph with vertices $\{1,\dots,n\}$ and edges $\{(i_k, j_k): 1 \leq k \leq d\}$ is a rooted forest. When $d = n-1$ this is the requirement that $\{(i_k, j_k): 1 \leq k \leq d\}$ define a rooted tree. There are famously $n^{n-1}$ such trees (and hence such monomials) by Cayley's formula (which tells us there are $n^{n-2}$ labeled trees on $\{1,2,\dots,n\}$, each of which can be rooted in $n$ different ways).
- The coefficient of $x^{k-1}$ in $d(xI - A)$ is in general a sum of $|\mathcal{F}_k|$ monomials in the variables $a_{ij}$ (each monomial having degree $n-k$). As just observed, it is well known that $|\mathcal{F}_1| = n^{n-1}$; I do not know of a simple citation for where the numbers $|\mathcal{F}_k|$ (equivalently, the number of rooted forests on $\{1,2,\dots,n\}$ involving $n-k$ edges) have been computed in print, but this MO post provides express (if somewhat ugly) formulas for the number of labeled forests with $k$ components on $\{1,\dots,n\}$ (that is, counting the objects in $\mathcal{F}_k$ up to differences in choices of roots).
- One can use computer software to expressly compute and check all of the formulas above, although they involve large numbers of terms even when $n$ is as small as $3$ or $4$, and the large amount of symmetry in the problem makes it difficult (at least for me) to visually process, let alone independently evaluate, the resulting output.