Probability of landing on the nth stair.

Let $p_{n,k}$ be the probability of hitting the $n^\text{th}$ stair after exactly $k$ steps, and let $p_n$ be the probability of hitting the $n^\text{th}$ stair at some point. Then $p_{n,k}$ is the coefficient of $x^n$ in the polynomial $(\frac13x+\frac13x^2+\frac13x^3)^k$. Therefore, $p_n=\sum_{k=0}^\infty p_{n,k}$ is the coefficient of $x^n$ in $$ \begin{align} \sum_{k=0}^\infty \left(\tfrac13x+\tfrac13x^2+\tfrac13x^3\right)^k &=\frac1{1-(\tfrac13x+\tfrac13x^2+\tfrac13x^3)} \\&=\frac{1/2}{1-x}+\frac{1/4}{1+x/(1+i\sqrt{2})}+\frac{1/4}{1+x/(1-i\sqrt{2})} \end{align} $$ The last equality can be derived using the usual partial fractions method. Expanding out each of these fractions into their Taylor series, you get $$ \boxed{p_n = 1/2+\frac{1/4}{(-1-i\sqrt{2})^n}+\frac{1/4}{(-1+i\sqrt{2})^n}} $$ So, you get an exact answer involving complex numbers magically canceling, along with an asymptotic result that the probabilities converge exponentially quickly to $1/2$.

Let $N(n,k)$ be the number of ways of reaching stair $n$ in exactly $k$ steps. If $P_n$ is the probability of reaching stair $n$, then:

$$P_n = \sum_{k = 0}^{\infty}3^{-k}N(n,k)$$

Note that $N(n,k)$ is defined for all integers $n$ and all nonnegative integers $k$, and so $P_n$ is defined for all integers $n$.

If you've reached stair $n$ in $k$ steps, and $k \geq 1$, that means you were at either stair $n-1, n-2$, or $n-3$ in $k-1$ steps. So:

$$N(n,k) = N(n-1,k-1) + N(n-2,k-1) + N(n-3,k-1); \qquad k \geq 1$$

Now, assume $n > 0$. We have:

\begin{align*} P_{n} &= N(n,0) + \sum_{k = 1}^{\infty}3^{-k}\big(N(n-1,k-1) + N(n-2,k-1) + N(n-3,k-1)\big) \\ &= \frac{1}{3}\sum_{k = 0}^{\infty}3^{-k}N(n-1,k) + 3^{-k}N(n-2,k) + 3^{-k}N(n-3,k) \\ &= \frac{1}{3}\big(P_{n-1} + P_{n-2} + P_{n-3}\big) \end{align*}

Since $N(n,0) = 0$ for $n \neq 0$. So, outside of $P_0 = 1$, each probability is just the average of the three prior probabilities, and $P_n = 0$ for $n < 0$. Some initial values:

$$P_0 = 1$$

$$P_1 = \frac{1}{3}(0 + 0 + 1) = \frac{1}{3}$$

$$P_2 = \frac{1}{3}(0 + 1 + 1/3) = \frac{4}{9}$$

$$P_3 = \frac{1}{3}(1 + 1/3 + 4/9) = \frac{16}{27}$$

$$P_4 = \frac{1}{3}(1/3 + 4/9 + 16/27) = \frac{37}{81}$$

Which agrees with your computation.

We have a Markov process

$\begin{bmatrix}p_{n+3}\\p_{n+2}\\p_{n+1} \end{bmatrix} =\begin{bmatrix}\frac {1}{3}&\frac 13&\frac 13\\1&0&0\\0&1&0 \end{bmatrix}\begin{bmatrix}p_{n+2}\\p_{n+1}\\p_{n} \end{bmatrix}$

$p_{n+2} = A^n p_0$

$A$ is diagonalizable, but I am getting a bit of a mess when I diagonlize it.

The eigenvalues are $1, \frac {-1+\sqrt 2 i }{3}, \frac {-1-\sqrt{2} i}{3}$