Let $P=\{p_{ij}\}$ be the transition matrix of an ergodic DTMC (say finite states for simplicity), let $p^{(k)}_{ij}$ be the entries of $P^k$.
If the associated Markov chain is initially at some state $i$ (i.e., $X_0=i$), and $V_n^{(j)}$ is the number of visits to some state $j$ (may not be $i$) up to time $n$, then it is true that $$ E[V_n^{j}|X_0=i]=\sum_{k=1}^np^{(k)}_{ij}. $$ This post has a related proof.
Also, is there any way to find the variance of $V_n^{j}|X_0=i$ using $P$ and (k-steps) transition matrices $P^k$, i.e, $$ Var[V_n^{j}|X_0=i] = ? $$
Thanks in advance!
$\newcommand{\1}{\mathbb{I}}$Let us call $\1$ the indicator function, i.e., $$\1_A(\omega) = \begin{cases}1, & \omega \in A\\ 0, & \omega \notin A.\end{cases}$$ Then we use two well-known facts, $$E(\1_A) = P(A), \, \text{Var}(\1_A) = P(A)(1-P(A)). \tag{1}$$ Now, we know that the number of visits is simply the sum of indicator functions, that is $$V_n^j = \sum_{k=0}^n\1_{X_k = j}. \tag{2}$$ Note, I took the sum from $0$ to $n$, but from the question, it is possible you want to take it to be $1$ to $n$. Now, taking expectations, we see that $$\begin{align} E(V_n^j) &= E\left(\sum_{k=0}^n\1_{X_k = j}\right) \tag{by (2)}\\ &= \sum_{k=0}^nE(\1_{X_k = j}) \tag{by linearity} \\ &= \sum_{k=0}^nP(X_k = j). \tag{by (1)} \end{align}$$ To get the variances, we just use the other identity in $(1)$.