Calculating probability of results of a non-deterministic function I am working on a computer program. It can be represented as a linear series of states $s_0 \ldots s_n$
Function $N(s_x)$ is defined such that:


*

*With a probability $p$ it will return the next state
$$s_x \rightarrow
\left\{
 \begin{array}{ll}
  s_{x+1} & \mbox{if } x \lt n \\
  s_{x} & \mbox{otherwise}
 \end{array}
\right.$$

*With a probability of $p$ it will return the previous state
$$s_x \rightarrow
\left\{
 \begin{array}{ll}
  s_{x-1} & \mbox{if } x \gt 0 \\
  s_{x} & \mbox{otherwise}
 \end{array}
\right.$$

*With a probability of $1-2p$ it will return the current state
$$s_x \rightarrow s_x$$


The program is an iterative application of $N$ to the initial value $s_i$ $r$ times $$Program = N^r(s_i)$$
I need to find $P(N^r(s_i) = s_x)$ for all $0 \leq x \leq n$ (the probability that applying $N$ iteratively to $s_i$ $n$ times yields the state $s_x$)
For small values I've been manually working out the different combinations that can lead to a particular state. For example, if $n=2$, $i=1$, $p=0.25$ and $r=2$ I can create a tree to help figure out the probabilities:

In this case:
$$P(s_0) = 0.25*0.25 + 0.25*0.5 + 0.5*0.25 = 0.3125$$
$$P(s_1) = 0.25*0.25 + 0.5*0.5 + 0.25*0.25 = 0.375$$
$$P(s_2) = 0.5*0.25 + 0.25*0.5 + 0.25*0.25 = 0.3125$$
However, this obviously gets impractical for the larger values I need to deal with in real life (e.g. $r=500$). Is there an easier way to calculate the probabilities?
 A: You've got a time-homogeneous Markov chain with finite state space, and you can diagonalize the transition matrix and decompose your initial state into its eigenvectors; then applying $N$ $r$ times just becomes multiplying each component with the $r$-th power of the corresponding eigenvalue. Of course this will only be more efficient than computing the $r$ steps directly if $n \ll r$.
A: Your probability is flowing according to the heat equation, and your boundary conditions are equivalent to a mirror, so if $i$=20 and $n$=99, then $s_x$ is the same as for an infinite line of heat flow initialized with point hotspots at ...,-221,-180,-21,20,179,220,279,... .
The heat at a point is just the sum of the contributions from these spots, and the contribution from a spot can be quickly estimated as $\frac{\exp{(\frac{-d^2}{4pr})}}{\sqrt{4\pi pr}}$, where $d$ is the distance from the heat source to the point of measurement.  These contributions die off very quickly (although any spot can contribute in principle if the time $r$ gets big enough), so you don't need to include spots from far away in your estimate.
