I suspect the confusion may be in what you mean by "step" in the phrase "$N$ steps ahead".
If you mean the one-second time step of your discrete time Markov process, then what you're calculating is exactly correct — with a $0.98$ self-loop transition probability, the process has a less than $0.5$ probability of moving anywhere within $N < 35$ steps.
If, however, by "step" you mean a transition in which the process does not stay in the same state, then you need to eliminate the self-loop transitions from your matrix (since you don't count them as steps) by setting their probability to $0$ and normalizing the remaining transition probabilities for each state so that they again sum to $1$.
Of course, what you lose by doing so is the information about the residence time in each state. If you wanted to know where the process will be after $N$ non-self-loop steps and how long it'll take to get there, that would require a rather more complicated calculation, and the resulting discrete phase-type distribution would, in general, be unlikely to have a simple form (although it might be possible to approximate it with some nicer distributions).
In the comments, you write that your transition matrix is (approximately) $$P=\pmatrix{0.99&0.0208&0&0\\0.0019&0.9792&0.0104&0\\0&0&0.9896&0.0108\\0.0028&0&0&0.9907}$$ and that "second option it is what I'm looking for", i.e. that you want to calculate the probability distribution after $N$ steps, a step denoting a transition between two distinct states. If so, what you need to do is to calculate another stochastic matrix $P'$ which describes the same process, but which does not include the self-loop transitions (i.e. whose diagonal entries are zero).
To obtain $P'$ from $P$, simply zero out the diagonal entries and re-normalize each row so that the matrix stays stochastic (i.e. so that each row sums to $1$). Since most of the rows in $P$ only have one non-zero off-diagonal element to begin with, that element becomes $1$, leaving row 2 as the only non-trivial case: $$P' = \pmatrix{0&1&0&0\\0.1545&0&0.8455&0\\0&0&0&1\\1&0&0&0}$$ From this, you can calculate the state distribution $\pi_n = \pi P'^n$ after $n$ steps, given a starting distribution $\pi$.
However, note that, whereas the original transition matrix $P$ was ergodic, $P'$ isn't — in particular, a process in an even-numbered state always ends up in an odd-numbered state after one step and vice versa, so all states have an even period. This lack of ergodicity doesn't really affect the calculation of $\pi_n$ as such, but it does mean e.g. that $\pi_n$ won't generally converge to a stationary distribution as $n \to \infty$.
Since $X$ is irreducible, for any state $j$ there exist positive integers $n,n'$ such that $P_{ij}^n>0$ and $P_{ji}^{n'}>0$. Since $P_{ii}>0$, it follows that $P_{ii}^m>0$ for all positive integers $m$, and hence $$P_{jj}^{n+n'+m}\geqslant P_{ji}^{n'}P_{ii}^mP_{ij}^n>0. $$ This implies that the period of state $j$ is $1$. Since $j$ was arbitrary, we conclude that $X$ is aperiodic.
Best Answer
This is my answer to my own question. I wanted to get an expression easy to work with, but this is the best I got.
The question is to compute the probability $P(N_1 \geq N | L)$ that is the probability that the number of visits to state $s_1$ (denoted by $N_1$) is larger than $N$ after $L$ state transitions. If state $s_1$ is revisited from state $s_1$ this is counted as an another visit. I count the visits to state $s_1$ as the number of state transitions that start from $s_1$; that is, the number of times that transitions $s_1 \rightarrow s_0$ and $s_1 \rightarrow s_0$ occur.
I denote by $P_1(n,l|s_i)$ the probability of visiting exactly $n$ times state $s_1$ in $l$ state transitions when the initial state is $s_i$. Then $P(N_1 \geq N | L)$ is \begin{equation} P(N_1 \geq N | L)=p_1\sum_{i=N_1}^{L}P_1(i,L|s_1)+p_0\sum_{i=N_1}^{L-1}P_1(i,L|s_0), \end{equation}
where the second term sums up to $L-1$ only because the initial state is $s_0$. I now compute $P_1(i,l|s_1)$
\begin{equation} P_1(i,l|s_1)=\sum_{j=1}^{l-i+1} \underbrace{p_{11}^{i-j}p_{10}^{j}\binom{i-1}{j-1}}_{A0} \underbrace{p_{01}^{j-1}p_{00}^{l-i-j+1}\binom{l-i}{j-1}}_{B0}+ \sum_{j=1_{l>i}}^{l-i} \underbrace{p_{11}^{i-j}p_{10}^{j}\binom{i}{j}}_{A1}\underbrace{p_{01}^{j}p_{00}^{l-i-j}\binom{l-i-1}{j-1}}_{B_1}, \end{equation}
where in $A0$-$B0$ and $A1$-$B1$ the final state is assumed to be $s_0$ and $s_1$, respectively. First I explain $A0$-$B0$. To visit state $s_1$ a total of $i$ times is it necessary to have exactly a total of $i$ $s_1 \rightarrow s_0$ and/or $s_1 \rightarrow s_1$ state transitions. In particular, there are $i-j$ and $j$ $s_1\rightarrow s_1$ and $s_1 \rightarrow s_0$ state transitions, respectively. Since the final state is $s_0$, and $s_1 \rightarrow s_0$ transition occurs $j$ times, there must be a total of $j-1$ $s_0 \rightarrow s_1$ state transitions. The total number of state transitions is $l$, then the total number of $s_0 \rightarrow s_0$ state transitions must be $l-j-i+1$. Now lets focus on $A0$, for each $s_0 \rightarrow s_1$ transition a set of $s_1\rightarrow s_1$ state transitions may occur (note 1), the total number of $s_1\rightarrow s_1$ transition sets is $j$ (i.e., number of $s_0 \rightarrow s_1$ transitions plus one). In total there must be $i-j$ $s_1\rightarrow s_1$ state transitions spread over $j$ $s_1\rightarrow s_1$ transition sets, the total number of possible combinations is the is the weak composition (note 2) of $j$ naturals that sum up to $i-j$ (note 3). $B0$ is obtained similarly, for each $s_1\rightarrow s_0$ transition there might be a set of $s_0\rightarrow s_0$ transitions. Hence there are a total $l-i-j+1$ $s_0\rightarrow s_0$ transitions spread over $j$ $s_0\rightarrow s_0$ transitions sets. Hence, the total number of possible combinations is the weak composition of $j$ natural numbers that sum up to $l-i-j+1$.
Now I will explain $A1$-$B1$. In this case the final state is $s_1$, and the number of $s_1 \rightarrow s_0$ state transitions is the same as the number of $s_0 \rightarrow s_1$ state transitions. The total number of $s_1 \rightarrow s_0$ and $s_0 \rightarrow s_1$ state transitions is $j$, the total number of $s_1 \rightarrow s_1$ is $i-j$, and the total number of $s_0 \rightarrow s_0$ is $l-j-i$. The total number of $s_1 \rightarrow s_1$ state transition sets is $j+1$, and the total number of $s_0 \rightarrow s_0$ state transition sets is $j$. The total number possible combinations of having $i-j$ $s_1 \rightarrow s_1$ state transitions spread over $j+1$ transitions sets is the weak combination of $j+1$ natural numbers that add to $i-j$. Following the same rationale, the total number possible combinations of having $l-j-i$ $s_0 \rightarrow s_0$ state transitions spread over $j$ transitions sets is the weak combination of $j$ natural numbers that add to $l-j-i$.
To compute $P_1(i,l|s_0)$ we can make use of $P_1(i,l|s_1)$,
\begin{equation} P_1(i,l|s_0)=\sum_{j=1}^{l-i}p_{00}^{j-1}p_{01}P_1(i,l-j|s_1), \end{equation}
where the term $p_{00}^{j-1}p_{01}$ corresponds to first transition probability from $s_0$ to $s_1$.
It is still needed to check all possible cases, and the limits of the sums. However, I think that this is pretty much it.
Notes